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PREFACE 


The  research  described  in  this  report  was  authorized  and  funded  under 
Work  Unit  32969,  “Development  of  a  New  Generation  Finite  Element  Harbor 
Wave  Model,”  of  the  Coastal  Research  Program,  sponsored  by  Headquarters, 

U.S.  Army  Corps  of  Engineers.  Administrative  responsibility  is  assigned  to  the 
U.S.  Army  Engineer  Waterways  Experiment  Station  (WES)  Coastal  and 
Hydraulics  Laboratory  (CHL),  Dr.  James  R.  Houston,  Director,  and  Mr.  Charles 
C.  Calhoun,  Jr.,  Assistant  Director.  Ms.  Carolyn  M.  Holmes,  CHL,  was  the 
Program  Manager,  and  Dr.  Zeki  Demirbilek,  Navigation  and  Harbors  Division 
(NHD),  CHL,  was  the  Principal  Investigator  for  the  work  unit.  Dr.  Demirbilek 
worked  under  the  administrative  supervision  of  Dr.  Martin  C.  Miller,  Chief, 
Coastal  Hydrodynamics  Branch,  and  Mr.  Gene  E.  Chatham,  Chief,  NHD,  CHL. 

This  report  describes  the  research  that  has  culminated  in  a  state-of-the-art, 
general-purpose  wave  predictive  model  called  CGWAVE.  Issues  that  are 
emphasized  in  this  report  include  theory  and  numerical  implementation 
aspects  of  this  model  and  a  set  of  examples  that  illustrate  the  application  of 
CGWAVE  to  real-world  problems.  A  step-by-step  user’s  guide  is  also  provided 
in  Section  7  to  facilitate  the  model’s  usage  in  projects. 

The  study  was  performed  and  the  report  prepared  over  the  period  1  June 
1995  through  15  August  1998.  Dr.  Demirbilek  and  a  team  of  researchers  led 
by  Dr.  Vijay  Panchang  from  the  University  of  Maine,  Orono,  Maine,  developed 
the  numerical  modeling  goals,  concepts,  and  methodology.  The  University  of 
Maine  team  included  Drs.  Bingyi  Xu  and  David  Stewart,  Messrs.  Liuzhi  Zhao, 
Karl  Schlenker,  Nishchey  Chhabra,  and  Wei  Chen.  The  study  team  completed 
development,  implementation,  and  testing  of  the  model.  The  field  validation 
part  of  this  research  was  collaborated  with  Dr.  Michele  Okihiro  and  Professor 
Robert  Guza  of  the  Scripps  Institution  of  Oceanography,  La  Jolla,  California. 

At  the  time  of  publication  of  this  report,  Director  of  WES  was  Dr.  Robert  W. 
Whalin.  Commander  was  COL  Robin  R.  Cababa,  EN. 


The  contents  of  this  report  are  not  to  be  used  for  advertising,  publication , 
or  promotional  purposes .  Citation  of  trade  names  does  not  constitute  an 
official  endorsement  or  approval  of  the  use  of  such  commercial  products. 
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Conversion  Factors, 

Non-SI  to  SI  Units  of  Measurement 


Non-SI  units  of  measurement  used  in  this  report  can  be  converted  to  SI  units  as  follows: 


Multiply 

By 

To  Obtain 

Acres 

4,046.873 

square  meters 

Degrees  (angle) 

0.01745329 

radian 

Feet 

0.3048 

meters 

Knots  (international) 

0.5144444 

meters  per  second 

Miles  (US  Statute) 

1.6093 

kilometers 

Nautical  miles 

1.852 

kilometers 

1  INTRODUCTION 


Wave  climate  plays  a  very  important  role  in  all  coastal  projects.  However,  in 
most  cases,  little  (if  any)  wave  data  are  available  for  engineering  construction  and 
planning.  Field  observation  and  physical  modeling  of  waves  are  extremely  difficult, 
costly,  and  time-consuming.  Buoys  are  far  away  from  the  project  site,  and  remote-sensing 
instruments  do  not  systematically  provide  wave  data  at  the  desired  resolution  in  the  near 
shore  region.  Since  no  data-recording  instrument  can  anticipate  future  sea  states,  the 
desired  sea-state  information  may  be  obtained  and  plans  evaluated  with  reliable 
mathematical  modeling  techniques. 

It  is  essential  to  have  reliable  information  on  wave  conditions  for  many  coastal 
and  ocean  engineering  problems.  The  most  important  wave  conditions  for  design  and 
assessment  in  project  studies  in  the  area  of  interest  include  the  wave  heights,  wave 
periods  and  the  dominant  wave  propagation  directions.  Typically,  these  wave  parameters 
are  obtained  from  a  wave  transformation  model  that  transfers  the  wave  data  collected  at 
some  remote  deep  water  site  to  the  location  of  the  project  in  the  near  shore.  As  waves 
move  from  deeper  waters  to  approach  the  shore,  these  fundamental  wave  parameters  will 
change  as  the  wave  speed  changes  and  wave  energy  is  redistributed  along  wave  crests  due 
to  the  depth  variation  between  the  transfer  sites  and  the  presence  of  islands,  background 
currents,  coastal  defense  structures,  and  irregularities  of  the  enclosing  shore  boundaries 
and  other  geological  features.  Waves  undergo  the  severest  change  inside  the  surf  zone 
where  wave  breaking  occurs  and  in  the  regions  where  reflected  waves  from  coastline  and 
structural  boundaries  interact  with  the  incident  waves. 

Until  recently,  the  linear  wave  ray  theory  was  used  for  wave  transformation  by 
tracing  rays  from  deep  water  to  the  project  site  near  shore.  The  effects  on  wave 
propagation  of  the  wave  height  and  direction  along  the  wave  crest  are  ignored  in  the  ray 
theory  since  this  theory  assumes  that  wave  energy  propagates  only  along  a  ray  and  thus. 
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energy  flux  is  conserved  between  two  adjacent  rays.  As  a  consequence  of  this 
assumption,  ray  theory  breaks  down  when  wave  ray  crossings  and  caustics  occur  because 
the  physics  of  diffraction  are  totally  ignored  in  the  numerical  ray  models. 

Starting  in  the  early  1980’s,  coastal  designers  and  researchers  have  recognized  the 
importance  of  the  combined  effects  of  refraction  and  diffraction  and  begun  to  develop 
improved  theories  and  associated  numerical  models.  There  are  indeed  several  wave 
theories  available  that  could  adequately  describe  the  combined  refraction  and  diffraction 
of  waves  from  deep  water  to  shallow  water  (Demirbilek  and  Webster  1992  and  1998). 
One  of  these  is  the  mild-slope  equation  (MSE).  This  is  a  depth-averaged,  elliptic  type 
partial  differential  equation  which  ignores  the  evanescent  modes  (locally  emanated 
waves)  and  assumes  that  the  rate  of  change  of  depth  and  current  within  a  wavelength  is 
small,  hence  the  ‘mild-slope’  acronym. 

Numerous  MSE-based  numerical  models  have  been  developed  for  predicting  the 
wave  forces  on  offshore  structures  and  studying  wave  fields  around  the  offshore  islands. 
Numerical,  laboratory  and  filed  tests  of  the  NSE  models  have  shown  that  the  MSE  can 
provide  accurate  solutions  to  problems  where  the  bottom  slope  is  up  to  1:3.  From  a 
practical  standpoint,  the  computational  requirements  for  solving  the  MSE  are  munch 
larger  than  those  for  ray  tracing.  The  reasons  for  this  are  because  the  MSE  is  a  two- 
dimensional  equation  and  has  to  be  solved  as  a  boundary-value  problem  with  appropriate 
boundary  conditions.  The  entire  domain  of  interest  must  be  discretized  and  solved 
simultaneously  and  the  element  size  has  to  be  small  enough  that  there  are  about  10  to  15 
nodes  within  each  wavelength.  These  requirements  place  severe  demands  on  computer 
resources  when  applying  MSE  models  to  large  coastal  domains. 

A  difficult  problem  in  the  prediction  of  waves  near  shore  is  to  determine  where 
approximately  the  wave  breaking  (and  breaker  line)  occurs  when  waves  are  inside  the 
sure  zone.  In  numerical  models  presently  used,  this  location  is  not  known  a  priori,  and  is 
usually  selected  with  an  ad  hoc  criteria  based  on  the  ratio  of  wave  height  to  local  water 
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depth.  Bottom  friction  and  dissipation  from  the  surrounding  land  boundaries  (i.e. 
entrance  losses  at  the  mouth  of  a  harbor)  may  also  be  empirically  incorporated  into  MSE 
models.  A  simplified  version  of  the  MSE  is  known  as  the  ‘parabolic  approximation’ 
(PA),  which  usually  greatly  reduces  the  excessive  computational  demands  of  MSE  model 
at  the  expense  of  further  assumptions  and  simplifications  which  may  render  the  numerical 
predictions  inaccurate  and  inappropriate  for  many  coastal  and  ocean  engineering 
problems  (Panchang  et  al.  1998). 

The  only  purpose  of  adapting  the  PA  is  to  convert  the  MSE  to  a  set  of  simpler 
equations  that  describe  a  wave  propagating  in  a  prescribed  direction  while  still  taking 
both  refraction  and  diffraction  in  the  lateral  direction  into  account.  The  greatest 
advantage  of  PA  is  its  numerical  efficiency,  it  can  be  solved  rather  easily  by  numerical 
means  and  thus  could  be  used  for  predicting  wave  transformation  over  a  relatively  large 
coastal  region.  When  reflection  is  of  major  interest,  as  it  is  in  harbors,  the  MSE  should 
be  used  since  the  PA  ignores  reflection.  One  must  also  be  reminded  that  the  PA  assumes 
that  the  length  scale  of  the  wave  amplitude  variation  in  the  direction  of  wave  propagation 
(x  direction)  is  much  longer  than  that  in  the  transverse  direction  (y  direction).  The  PA  is 
derived  on  the  assumption  that  percentage  changes  of  depth  within  a  typical  wavelength 
are  small  compared  to  the  wave  slope.  For  details  about  PA  models,  see  Booij  (1981), 
Liu  (1983),  Kirby  (1983),  Liu  and  Tsay  (1984),  and  Kirby  and  Dalrymple  (1984).  The 
PA  has  been  verified  extensively  by  laboratory  studies  and  field  applications  (Berkhoff  et 
al.  1982),  Liu  and  Tsay  (1984),  Kirby  and  Dalrymple  (1984),  Vincent  and  Briggs  (1989), 
(Demirbilek  1994,  Demirbilek  et  al.  1996a  and  1996b),  and  Panchang  et  al  (1998). 

The  mild-slope  wave  equation  (also  known  as  the  "combined  refraction- 
diffraction"  equation),  first  suggested  by  Eckart  (1952)  and  later  re-derived  by  Berkhoff 
(1972,  1976)  and  others,  is  now  well-accepted  as  the  method  for  estimating  coastal  wave 
conditions.  It  can  be  used  to  model  a  wide  spectrum  of  waves,  since  it  passes,  in  the 
limit,  to  the  deep  and  shallow  water  equations.  Although  the  equation  was  developed  in 
the  mid-seventies,  computational  difficulties  precluded  the  development  of  a  model  for 
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the  complete  mild-slope  equation  (except  for  very  small  domains).  Typically,  coastal 
wave  propagation  problems  involve  the  modeling  of  very  large  domains.  For  example, 
consider  the  case  of  12  second  waves  in  water  of  15  m  depth.  The  wavelength  L  is  about 
136  m;  an  8  km  by  8  km  domain  is  about  3600L2  in  size.  The  difficulties  associated  with 
solving  such  large  problems  spawned  the  development  of  several  simplified  models  (e.g. 
the  "parabolic  approximation"  models  (Dalrymple  et  al.  1984;  Kirby,  1986),  RCPWAVE 
model  (Ebersole,  1985),  EVP  model  (Panchang  et  al  1988),  etc.).  However,  these 
simplified  models  compromised  the  physics  of  the  mild-slope  equation:  they  model  only 
one-  or  two-way  propagation  with  weak  lateral  scattering.  Such  models  are  hence 
applicable  only  to  rectangular  water  domains  for  a  very  limited  range  of  wave  directions 
and  frequencies.  Most  realistic  coastal  domains  with  arbitrary  wave  scattering  cannot  be 
modeled  with  these  simplified  models. 

This  manual  describes  a  wave  model  called  CGWAVE  developed  at  the 
University  of  Maine  under  a  contract  for  the  U.S.  Army  Corps  of  Engineers,  Waterways 
Experiment  Station.  CGWAVE  is  a  general  purpose,  state-of-the-art  wave  prediction 
model.  It  is  applicable  to  estimation  of  wave  fields  in  harbors,  open  coastal  regions, 
coastal  inlets,  around  islands,  and  around  fixed  or  floating  structures.  While  CGWAVE 
simulates  the  combined  effects  of  wave  refraction-diffraction  included  in  the  basic  mild- 
slope  equation,  it  also  includes  the  effects  of  wave  dissipation  by  friction,  breaking, 
nonlinear  amplitude  dispersion,  and  harbor  entrance  losses.  CGWAVE  is  a  finite-element 
model  that  is  interfaced  to  the  SMS  model  (Jones  &  Richards,  1992)  for  graphics  and 
efficient  implementation  (pre-processing  and  post-processing).  The  classical  super¬ 
element  method  as  well  as  a  new  parabolic  approximation  method  developed  recently 
(Xu,  Panchang  and  Demirbilek  1996),  are  used  to  treat  the  open  boundary  condition.  An 
iterative  procedure  (conjugate  gradient  method)  introduced  by  Panchang  et  al  (1991)  and 
modifications  suggested  by  Li  (1994)  are  used  to  solve  the  discretized  equations,  thus 
enabling  the  modeler  to  deal  with  large  domain  problems.  This  manual  provides  a  brief 
review  of  the  basic  theory  in  Sections  2  and  3,  an  overview  of  how  this  theory  is 


4 


implemented  in  Sections  4  through  6,  a  step  by  step  guide  for  using  this  wave  model  in 
Section  7,  and  a  set  of  examples  that  were  run  using  CGWAVE  in  Section  8. 
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2  BASIC  EQUATIONS 


The  solution  of  the  two-dimensional  elliptic  mild-slope  wave  equation  is  a  well- 
accepted  method  for  modeling  surface  gravity  waves  in  coastal  areas  (e.g.  Chen  & 
Houston,  1987;  Chen,  1990;  Xu  &  Panchang,  1993;  Mei,  1983;  Berkhoff,  1976;  Kostense 
et  al.,  1986;  Tsay  and  Liu,  1983).  This  equation  may  be  written  as: 


where 
fj(x,y)  = 


a  = 

C(x,y)  = 
Cg(x,y)  = 


k(x,y)  = 


complex  surface  elevation  function,  from  which  the  wave 
height  can  be  estimated 

wave  frequency  under  consideration  (in  radians/second) 

phase  velocity  =  cr/k 

group  velocity  =  do  /  8k  =nC  with 

1 r  2kd  ^ 

n  =  —  1  + - 

2  ^  sinh  2kd  y 

wave  number  (=  2it/L),  related  to  the  local  depth  d(x,y) 
through  the  linear  dispersion  relation: 

=  gk  tanh  (kd) 


(1) 


(2) 


(3) 


Equation  1  simulates  wave  refraction,  diffraction,  and  reflection  (i.e.  the  general 
wave  scattering  problem)  in  coastal  domains  of  arbitrary  shape.  However,  various  other 
mechanisms  also  influence  the  behavior  of  waves  in  a  coastal  area.  The  mild-slope 
equation  can  be  modified  as  follows  to  include  the  effects  of  frictional  dissipation 
(Dalrymple  et  al  1984;  Chen  1986;  Liu  and  Tsay  1985)  and  wave  breaking  (Dally  et  al 
1985;  De  Girolamo  et  al  1988): 
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(4) 


+  iaw  +  iC  ay 


fj  =  0 


where  w  is  a  friction  factor  and  y  is  a  wave  breaking  parameter.  Following  Dalrymple  et 
al-  (1984),  we  have  used  the  following  form  of  the  damping  factor  in  CGWAVE: 


2ntf>| 

2fr  ak2 

k  J 

37t  (2kd  +  sinh  2kd)  sinh  kd 

(5) 


where  a  (=  H/2)  is  the  wave  amplitude  and  fr  is  a  friction  coefficient  to  be  provided  by  the 
user.  The  coefficient  fr  depends  on  the  Reynolds  number  and  the  bottom  roughness  and 
may  be  obtained  from  Madsen  (1976)  and  Dalrymple  et  al.  (1984).  Typically,  values  for 
fr  are  in  the  same  range  as  for  Manning’s  dissipation  coefficient  ‘n\  Specifying  fr  as  a 
function  of  (x,y)  allows  the  modeler  to  assign  larger  values  for  elements  near  harbor 
entrances  to  simulate  entrance  loss.  For  the  wave  breaking  parameter  y,  we  use  the 
following  formulation  (Dally  et  al  1985,  Demirbilek  1994,  Demirbilek  et  al.  1996b): 


y 


xfi  r2d2) 
d l  4a2  J 


(6) 


where  %  is  a  constant  (a  value  of  0.15  is  used  in  CGWAVE  following  Dally  et  al  (1985)) 
and  T  is  an  empirical  constant  (a  value  of  0.4  is  used  in  CGWAVE). 

In  addition  to  the  above  mechanisms,  nonlinear  waves  may  be  simulated  in  the 
MSE.  This  is  accomplished  by  incorporating  amplitude-dependent  wave  dispersion, 
which  has  been  shown  to  be  important  in  certain  situations  (Kirby  and  Dalrymple  1986). 
The  nonlinear  dispersion  relation  used  in  place  of  Equation  3  is 
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<52  =  gk[l  +  (ka)2F,  tanh5  kd]tanh{kd  +  kaF2 } 


(7) 


where 


Fi  = 


f2  = 


cosh(4kd)  -  2  tanh2  (kd) 
8sinh4(kd) 

'  kd  y 

k  sinh(kd)  J 


(8) 
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3  BOUNDARY  CONDITIONS 


Along  rigid,  impermeable  vertical  walls,  no  flow  normal  to  the  surface  gives 
3f|  /  5n  =  0 .  However,  in  general,  the  following  partial  reflection  boundary  condition 
applies  along  coastlines  or  permeable  structures 


d  f| 

d  n 


=  ar| 


(9) 


where  a  =  a,  +ia2  is  a  complex  coefficient .  For  simplicity,  a  is  generally  represented 
as 


a  =  ik 


1-Kr 
1  +  Kr 


(10) 


where  Kr  is  the  reflection  coefficient  (Tsay  and  Liu,  1983;  Chen  and  Houston  1987). 

Along  the  open  boundary  where  outgoing  waves  must  propagate  to  infinity,  the 
Sommerfeld  radiation  condition  applies 


lim 
kr— >oo 


Mr 


ik 


Os  0 


(ID 


where  fis  is  the  scattering  wave  potential.  It  is  shown  in  Mei  (1983)  that  the  desired 
scattered  wave  potential  fjs ,  which  is  a  solution  of  the  mild-slope  equation  and  satisfies 
the  radiation  condition  Equation  11,  can  be  written  as: 


Os  =  2 cos  110  +  Pn  sin  n0) 


n=0 


(12) 
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where  H„(kr)  are  the  Hankel  functions  of  the  first  kind.  The  Hankel  functions  of  the 
second  kind  do  not  satisfy  the  Sommerfeld  radiation  condition  at  infinity  and  are  hence 
excluded  from  (12). 

However,  the  T|s  given  in  (12)  requires  that  the  exterior  domain  be  of  constant 
depth.  Also  for  harbor  problems  (Figure  1),  the  scattered  wave  potential  as  described  by 
(12)  demands  straight,  collinear  and  fully  reflective  coastlines  in  the  exterior  region.  To 
overcome  these  problems,  Xu,  Panchang  and  Demirbilek  (1996)  have  developed  an 
alternative  scheme  in  dealing  with  the  open  boundary  condition.  This  consists  of  using 
the  following  parabolic  approximation  along  the  open  boundary  : 


df)s 

— 12-  +  pris+q 
dr 


=  0 


(13) 


where 


k2r2  +k;jr2  +ik0r  +  - 


P  = 


2ik0r 


and  q 


2ik0r 


(14) 


In  Equation  14,  ko  can  be  taken  as  the  wave  number  corresponding  to  the  averaged  water 
depth  along  the  open  boundary  T.  Within  the  model  domain  Q.,  the  mild-slope  equation 
applies.  The  parabolic  approximation  (13)  will  be  used  only  along  the  semi-circular  arc  T 
as  the  open  boundary  condition.  The  actual  implementation  of  these  boundary  conditions 
is  described  later. 
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4  FINITE  ELEMENT  FORMULATION 


4.1  Open  Sea  Problems 

CGWAVE  uses  the  finite-element  method,  which  is  a  powerful  approach  for 
modeling  coastal  phenomena  in  regions  of  complex  shape.  In  the  case  of  a  group  of 
scatterers  surrounded  by  an  open  sea  of  constant  depth,  the  incident  wave  may  be  written 
as  (Demirbilek  and  Gaston,  1985) 


fjj  =  Aeikrcos(e-e,)  =  n(e  “  ei) 

n=0 


(15) 


where  A  is  the  amplitude  of  the  incident  wave,  6i  is  the  incident  wave  angle  with  respect 
to  the  x-axis,  Jn  is  the  n-th  order  Bessel  functions  of  the  first  kind  and 

fl  when  n  =  0 

8  =\  (16) 

“  [2  when  n  *  0 

The  incident  direction  is  defined  such  that  the  incident  wave  travels  in  the  positive  x- 
direction  when  0i  is  equal  to  zero;  the  x-direction  is  obtained  from  the  bathymetry  data 
file  that  is  input  to  the  CGWAVE  program.  This  bathymetry  file  should  be  oriented  such 
that  the  x-axis  points  to  the  east. 
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As  shown  in  Figure  2,  the  entire  domain  is 
separated  into  two  sub-domains.  Domain  Q  is  the 
numerical  model  domain.  Domain  Q0  is  the  exterior 
domain  extending  to  infinity.  We  assume  that 
complicated  topography,  structures,  and  islands,  are 
located  inside  the  circular  boundary  F  (in  domain  Q). 
In  Qo>  the  total  wave  potential  can  be  written  as  the 
sum  of  incident  wave  potential  and  the  scattered 
wave  potential: 


Fig.  2  Definition  sketch 


flex.  =r\i+f\s 


(17) 


For  brevity,  we  write  the  governing  Equation  4  in  the  general  form: 


V-(aVf|)+  bfj  =  0 


(18) 


where  a  =  CCg  and 


b  a 


~g2  +  iow  +  iCgOY  • 


Mei  (1983)  has  shown  that  the  problem  of  solving  Equation  18  with  boundary 
conditions  described  by  (9)  on  coastlines/structures  and  by  Equation  11  at  infinity  is 
equivalent  to  the  stationary  of  the  following  functional  J: 


J  = 


JJ|[a(vnf 


aafi2ds  + 


dn 


dn 


Ids 


(19) 


The  solution  of  the  wave  potential  can  be  found  by  minimizing  J  over  domain  Q. 
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In  the  finite-element  method,  we  first  discretize  the  computational  domain  Q  into 
a  network  of  simple  triangular  elements.  The  size  of  these  elements  should  be  much 
smaller  than  both  the  local  wavelength  and  the  scale  of  local  bathymetric  variation.  Finer 
resolution  is  also  desirable  at  places  where  the  change  of  amplitude  in  space  is  rapid  (e.g. 
near  “caustics”).  Over  each  triangular  element,  the  wave  potential  f\  is  approximated  by 

the  following  linear  two-dimensional  function  f) e , 


+ Nj + m 


(20) 


where  f|  ie  represent  the  wave  potentials  at  the  comers  (nodes)  of  the  element  e  (Figure  3) 
and  N®  are  the  linear  interpolation  functions: 


Ne  _  a,  +bix  +  ciy 
1  2Ae 


with 


a;  =xjyk-yjxk 

bi=yj-yk 

ci=xk-xj 


(22) 


and 


(21) 


Fig.  3  A  typical  element 


Ae  =  area  of  element  e 
*;  y'x 
=  21  x2  y2e 

1  x3e  y3 


(23) 
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Note  that  for  future  use,  we  give 
JjN^dxdy  = 


— Ae,  i  =  j  =  k 

60  J 
A 

fff  N*N-N£  dx  dy  =  \—Ae,  i  =  j  or  i  =  k  or  j=  k 

JJJ  60 

■j~Ae,  i*  j*k 


Notice  that  in  the  above  formulations,  (i,  j,  k)  are  denoted  in  a  counter-clockwise  manner. 


For  element  e,  the  following  relations  can  be  established  for  substitution  into  (19): 


Y7-'e  '•"V'1  '■  e  dN*  "V’  *■  e  dN* 

VTie  =  iX^  ^+jS7ii-r±- 

j— i  ox  j_j  dy 

(Vfle)2  =  Vfle-Vfie 


dNf  dN° 
dx,  dx. 


=/^iT  MM 

*•  dx,  dx. 


dNj  dNM  [  dN®  dN 


dxa  dxa  J  ^  dxa  dxa 


dN^  dNj 
dx,  dx. 


f  ?  M 
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(26) 


(Nr)2  (NfN;) 
(N^r)  (n^)2 

(N^Nf)  (NOT) 


where  a  =  1,  2  is  the  dummy-index  notation  and  (x,  =  x,  x2  =  y),  i.e.  Equation  25 
represents  the  sum  for  a  =  1  and  a  =  2.  Note  that  { f|e}T  =  [f|  *f\zsf\3e]  from  Equation 
20. 


We  may  also  assume  that  the  coefficients  (a,  b)  in  Equation  19  vary  linearly  on 
element  e: 


(27) 


For  the  first  part  of  Equation  19,  we  may  write 


n  ~ 

z  eeW  ]x3  3x3  3x1 


(28) 


where 


Ku.j 


( 


Ne 

1NP 


3N?  9N‘ 


dx„  8x 


dxdy  -  bpjjNpNfNj  dxdy 


<*  ) 


e 


(29) 
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where  (3  =  1, 2,  3  is  another  dummy-index  notation.  Since 


3Nf  3N- 
9xa  dxa 


(30) 


and 


II 


Npdxdy : 


we  have 


r  r 


NS 


dxa  dxa  f 


dxdy  = 


a,e  +  a,e  +  a,e 


12Ae 


(bibj+c^) 


(31) 


(32) 


The  second  term  of  Equation  29  is 


bpJjNpNfN-dxdy 

e 

=  bf  |f  (Nf  )2  N  ®  dxdy  +  bj  |f  Nf  (N  *  )2  dxdy  +  bj;  JJ  NfN  *  N£dxdy  (33) 

e  e  e 

=  —  (2bie+2bt+b") 

60 V  1  J  k/ 


for  i  *  j ,  and 


(34) 


SfjjNJNfNJdxdy 

e 

=  b'  j|(N*)3 dxdy  +  b*  }Jn;,(n')2 dxdy  +  bk'2 § NJ, (N* )2 dxdy 

e  e  e 

=  ^(3b'  +  b;,+bl%) 

for  i  =  j  =  i ,  where  kl  &  k2  are  the  other  two  nodes  of  element  e.  Now,  Equation  29  can 
be  written  as 


Ku.j  = 


a,  +  a2  +  a3 
12Ae 


(b,bj 


— (2b'+2b'+bj) 
60 V  1  J 


A l 
.30 


+  CiCj)- 
when  i  *  j 


when  i  =  j  =  i 


(35) 


After  computing  the  element  matrix  [K®]  for  all  the  elements  (e  =  1,  2,  3, ... ,  E),  where  E 

is  the  total  number  of  elements,  we  can  assemble  them  into  a  “global”  matrix  [Kj], 
Equation  28  becomes 

i.  =^Z{n'}T[Kf]{f«"}=|-{n}T[K,Kfi}  <36> 

et=n  1x3  3x3  3x1  ^ 


In  (36) 


{t\}T  =  {^1^  . ’  t)n} 


where  N  is  the  total  number  of  nodes  in  domain  Q 


(37) 
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The  second  part  of  Equation  19  is 


I2  =  f  ^-aafj2ds  (38) 

J  B  £ 

where  B  denotes  all  the  coastline  and  internal  land 
boundaries. 

Along  a  segment  p  of  the  coastal  boundaries  B 
(Figure  4),  the  wave  surface  elevation  function  fj  and 
the  variable  a  are  approximated  by  linear  functions  as 
before: 


ff  =  Nffjf  +  N?fjj  (39) 

Figure  4.  Boundary  segments 

ap  =N?ajp+NJpap  (40) 

Here,  we  assume  that  i  -»  j  is  the  positive  direction  of  the  boundary,  counter-clockwise 
for  the  coastline  boundaries  and  clockwise  for  internal  land  boundaries. 

First,  let  us  find  the  expression  for  ( fj p)2  on  segment  P: 


|T  K)2 

(NfN;)' 

IKnd 

(nO2  . 

AP 

J 

(41) 


Equation  38  can  then  be  written  as 


(42) 


I2  =  f  ^aatfds 

J  B  2* 


fV:1,  .  _/ 

-  P 

/ 


4£ftr.<tfK 

*  P=1 


where  Nb  is  the  total  number  of  nodes  along  boundary  B  and 


K^^XNfN'ds 

-i-cxLp(ap  +  ap)  for  i*j 
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— aLp (3ajP  +  ap)  for  i  =  j  =  i 

112  v  J  ' 


(43) 


Assembling  all  segments  on  coastal  boundary  B,  we  have 


Nh-1 


ftn  i 


i,=|l{Cn[}[K']  4ft'}  [Kjlft'} 

Z  lxNb  NbxNb  Nbxl 


v^j  J 


(44) 


The  third  part  integral  in  Equation  19  is 


dOs  fid(4+  Tji) 
3n  9n 

ds 


ds  = 


i3 + 14 + 15 + 16 


-J. 

=+J, 


ai)— — ds 
g  dn 

arj— ^ds 
g  dn 

aflj^-ds 
g  dn 


i5 

L 


(45) 
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For  simplicity,  the  open  boundary  T  is  assumed  to  be  a  circle  of  radius  R.  For 
computational  purposes,  the  series  for  the  scattered  waves  (Equation  12)  and  the  series 
for  incident  waves  (Equation  15)  are  truncated  after  a  finite  number  of  terms.  In 
principle,  trial  and  error  should  be  performed  in  modeling  a  certain  case  in  order  to 
choose  the  appropriate  number  of  terms.  Here,  we  assume  that  the  series  will  be  truncated 
after  m  terms. 


By  using  the  orthogonality  of  trigonometric  functions: 


n0  sin  m0d0 


0  when  n^m 
K  when  n  =  m 


(46) 


r  2n 


cosn0cosm0d0  = 


0 

K 

2n 


when  n^m 
when  n  =  m  ^  0 
when  n  =  m  =  0 


(47) 


and  substituting  f)  s  with  Equation  12,  the  line  integral  I3  in  Equation  45  can  be  evaluated 
analytically,  as  follows 


_ .  71 ,  ^  _ 

I,  =  — kRa 


m 

zocJhX  +£(°^  +PJ.)H„Hd 


n=l 


(48) 


where  k  and  a  can  be  taken  as  average  values  along  T  and 


H.-H®(kR),  Hn  = 


d(kr) 


H®(kr) 


(49) 


r=R 


For  convenience  of  mathematical  manipulation  later  on,  we  define  the  following  vector 
for  the  unknown  coefficients  a*  and  Pi : 
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(50) 


M  —  {^0’  ^1  ’  Pl»  ^2’  P29  ***  ’  ^m’  Pm} 

IxM 


where  M  =  2m  +1.  The  integral  I3  can  now  be  rewritten  as 


13  =^-{h}t[k!Kh} 


where  [K3]  is  a  diagonal  matrix  of  dimension  M  by  M: 


(51) 


[K3]  =  7tkRa  diag{2H0H0,  H,H1;  H,!!,,  ...  ,  HmHm,  HmHm}  (52) 


The  integral  I4  in  Equation  45  is 


-  dns 


I‘=LSfllfdS 

Nr 


P=1 

Nr 


=  ^1  (Nfftf  +  Nff,') 


segment  P 


ill 

a0H0  +  ZH  n(ancosn0p  +pDsinn0p)  Us 


n=0 


L  P=i 


111 

a0H0  +  lH.(a  ncosn0p  +Pnsinn0p) 


n=0 


(53) 


where  Lp  is  the  length  of  segment  P  and  Nr  is  the  total  number  of  segments  (=  total 
number  of  nodes)  along  the  circular  boundary  T  (Figure  5).  In  Equation  53,  the  value  of 
3  -q  s/3n  is  approximated  by  its  value  at  the  center  of  segment  P.  Equation  53  may  be 
written  in  matrix  form  : 


i4=fer7[K4M 


(54) 
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where  { fj r}  is  the  subset  of  { fj }  for 
nodes  situated  on  boundary  T : 

{nr}T  =  {nl\  C  -  .  r\lr} 

lxNr 

(55) 

and  [K4]  is  a  fully  populated 
Nr  xM  matrix  : 

[K,]  =  f  X 

2H0L1  •••  Hn(cosn0Nr_,  H-cosnB^L1 
2H0L2  •••  Hn(cosn0, +cosn02)L2 


pH0LNp  •  •  •  Hn (cosn0Nr_j  +  cosn0Nr  )ln’ 
where  n  =  1, 2, m. 


Figure  5.  Definition  sketch  for  line  integrals 
along  the  open  boundary  T 


H j  (sin  n0Nr_j  +  sin  n0,  )l‘  •  •  • 
Hn  (sin  n0,  +  sin  n02  )L2 


Hn  (sin  n0Nr_,  +  sin  n0Np  )l2 


(56) 

The  next  integral  in  Equation  45  is  I5.  Similar  to  the  treatment  in  I4,  we  will  take 
the  center  value  of  0^  i/3n  and  assume  linear  variation  of  fj  for  a  segment  on  boundary  T. 
Substituting  fj  1  in  I5  by  Equation  15,  it  is  easy  to  find 


1  Nr 

=  —  ka  A]£  Lp  (fif  +  r\p  )cos(0p  -  0T  )exp[ikRcos(0p  -  0, )]  (57) 

2  p=i 

={QjT{flr} 
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where 


+  q,)L\  (q,+q2)L2,  ...  ,  (qNr_: +qNr )LNr } 


(58) 


and 


qp  =  cos(0p  -0I)exp[ikRcos(0P  -0j)],  P  =  1, 2, ... ,  Nr  (59) 


The  last  integral  1«  in  Equation  45  involves  both  f|  s  and  fj  i.  Here,  the  Bessel- 
Fourier  form  of  f|i  (see  Equation  15)  will  be  used  and  the  integral  can  be  found 
analytically: 


=  kaA  (*  ^  £nin J n  (kr)cosn(0  -  0j  )£  Hn  (ancosn0  +  (3nsinn0)ds 

=0  n=0 

£eninJn  (kr)cosn(0  -  0,  )£  Hn  (ancosn0  +  pnsinn0)Rd0 


r  u=o 

m 


=  kaA 


(60) 


0  n=0 


={Q.m 


where 


{Q6}t  =  27tkRaA  x 

{j0H0,  iJjHicos©,,  iJ,H;sin0j,  ...  ,  imJmHmcosrn0I,  imJmHmsinrn0I} 


(61) 


Now,  we  have  evaluated  all  the  integrals  of  the  functional  J  defined  by  Equation 
19  using  a  linear  triangular  element  network.  Collecting  these  integrals  together,  we  have 
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(62) 


J  =  I,-I2+l3-I4-I5+I6 

4Wt[k,  ]{f,}  - 1  {f|- }T[K!]{f,' } + 1  >}t[k3  Mu} 

-{r}T[K4]{(l}-{Q,}T{fir}+{Q6}T{tl} 


Since  J  is  stationary,  the  following  must  be  true  for  the  solution  of  the  problem  : 


^  =  0  i  =  1, 2, 3, N  (63) 

drji 


and 


aj 

apj 


=o 


j=  1,2,3, 


M 


(64) 


These  relations  give 


(65) 


and 


[K3M-[K4]T{f)r}  =  -{Q6} 


(66) 


From  Equation  66,  we  have 


{m}=-[K3r{Q6}+[K3r[K4]T{fir} 


(67) 


25 


Substituting  {ji}  in  Equation  65  by  Equation  67,  Equation  65  becomes 


[K,]{f|}-[K,]{fi*}-[K4][K3nK4]T{fir}  =  {Qs}-[K4lK,r,{Q6}  (68) 

Finally,  after  proper  assembling,  we  have 


[A]{<)}  =  {f} 


(69) 


Equation  69  is  the  desired  linear  system  of  equations,  which  is  the  finite-element 
representation  of  the  mild-slope  equation  for  open  sea  problems.  Notice  that  the  boundary 
conditions,  including  coastline  boundaries  and  the  circular  open  boundary,  are  all 
consolidated  in  Equation  69.  The  solution  method  used  in  CGWAVE  for  Equation  69  is 
described  in  the  next  section. 

4.2  Harbor  Problems 

The  finite-element  formulation  given  above  is  for  open-sea  offshore  problems.  In 
case  of  harbor  problems,  the  formulation  is  analogous.  The  only  difference  arises  from 
the  treatment  of  the  open  boundary  condition.  The  classical  treatment  of  these  problems 
assumes  that  the  coastlines  outside  the  model  domain  are  straight,  collinear  and  fully 
reflective.  The  exterior  wave  field  is  written  as  f|  ext  =  fl  i  +  R  +  s>  where  f| , ,  fj  R,  and 
f|  s  represent  the  incident,  the  reflected,  and  the  scattered  wave  fields,  respectively.  Based 
on  the  assumptions,  we  define  (Demirbilek  and  Gaston  1985) 

^0=^1+ “Hr 

_  Aeikrcos(e-0, )  +  Aeikrcos(0+6, )  (70) 

oo 

=  2A^  eniD  Jn  (kr)cosn0jcosn0 

n=0 


where  A  is  the  incident  wave  amplitude  and  0i  is  the  incident  wave  angle  with  respect  to 
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the  exterior  coastlines  as  shown  in  Figure  1.  The  scattered  wave  potential  fjs  in  the 
exterior  region  must  take  the  following  form  in  order  to  comply  with  the  exterior 
coastline  boundary  conditions: 


f|s  =  ]£Hn(kr)ancosne 

n=0 


(71) 


as  shown  in  Xu,  Panchang  and  Demirbilek  (1995).  The  corresponding  functional  for 
harbor  problems  has  the  same  form  as  Equation  19  except  that  f)  i  in  Equation  19  has  to 
be  replaced  by  fjo  (Equation  70),  fjs  takes  the  new  form  given  by  Equation  71,  and  the 
open  boundary  T  represents  the  semicircle  as  shown  in  Figure  1.  The  finite-element 
formulation  of  harbor  problems  can  now  readily  be  found  in  a  manner  similar  to  the  open- 
sea  problems  described  above,  by  replacing  f\  \  and  fj  s  with  Equation  70  and  Equation  71 
and  performing  the  boundary  integration  for  I4  through  1$  from  0  to  K. 


4.3  Alternative  Open  Boundary  Treatment 

For  most  practical  cases,  the  fully  reflective  straight  coastline  assumption  in  the 
classical  treatment  of  the  open  boundary  condition  is  improper  and  the  effects  may 
substantial.  Xu,  Panchang  and  Demirbilek  (1995)  have  shown  that  it  is  preferable  to  use 
the  parabolic  approximation  (Equation  13)  as  the  open  boundary  condition.  For  harbor 
problems,  along  the  open  boundary  T  (Figure  1)  we  use  (13)  as  the  boundary  condition 
for  the  scattered  waves.  Matching  the  potential  and  its  normal  derivative  along  T  and 
using  the  parabolic  open  boundary  condition  (13),  we  have  the  total  potential  as 


and 


r\  =  f\0+  f\s 


afi  _  gfjo  [  afjs  _  af\0 

9n  9r  9r  9r 


P^ls+q 


aV| 

ae2 , 


(72) 


(73) 
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Using  Equation  72  to  eliminate  fj  s  yields 


an 

dn 


„  d2f\ 

'm~qW  g 


(74) 


where 


g  = 


ago 

3r 


+  Pfio+<J 


a2n0 

ae2 


(75) 


Equation  74  is  the  open  boundary  condition  in  terms  of  total  wave  potential  f\ . 

For  this  parabolic  boundary  problem,  the  desired  Jacobian  functional  may  be 
more  complicated  than  that  of  Equation  19.  Therefore  the  Galerkin  finite-element 
formulation  is  used. 


According  to  the  Galerkin  approach,  the  wave  potential  f)  is  approximated  by 


(76) 


where  f[  j  is  the  solution  fj  at  node  i  and  Nj  (x,  y)  is  the  linear  interpolation  function  for 
node  i.  The  unknown  f) ;  can  be  determined  from  the  orthogonality  conditions  between 
function  N;  and  left-hand  side  of  Equation  1 8,  that  is 

JJ  (V(aVq)+  bq^NjdQ  =  0  (i  =  1,...,N)  (77) 

Cl 

Using  the  divergence  theorem.  Equation  77  becomes 
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(78) 


JJ  aVfiVNjdQ  +  JJ  bfjVNjdn  =  0 

n  £2 


where  C  and  T  denotes  coastal  boundaries  and  open  boundaries  respectively,  and  dfj/dn 
is  the  normal  derivative  of  fj . 


Under  the  linear  assumption,  the  interpolation  function  N;  is  to  satisfy 


Nj(x,  y)  = 


'1 

io  to  1 


[0 


for  (x,  y)  at  node  i 

for  (x,  y)  within  elements  surrounding  node  i 
for  (x,  y)  outside  elements  surrounding  node  i 


(79) 


Therefore,  Ni  can  be  represented  as 


N,(x,y)  =  £NT(x,y)  (80) 

ei 

where  e*  refers  to  those  elements  around  node  i  and  N;e(x,y)  is  the  linear  interpolation 
function  corresponding  to  an  element  e  and  one  of  its  node  i.  When  (x,  y)  is  at  boundary, 
Equation  80  becomes 

Ni(x,  y)  =  £Nf(x,y)=N?(x,y)+Ni*(x,y)  (81) 

where  Pi  and  P2  are  the  boundary  segments  to  either  side  on  node  i.  The  function  Np=0 
for  all  other  segments. 

Within  an  element  e,  the  value  of  f|  is  obtained  by  substituting  (80)  into  (76), 
t)(x,  y)  =  'njNj  +  f|kNek  (82) 


29 


where  i(xj,y9,  j(Xj,yj),  k(xk,yk)  are  the  three  nodes  of  element  e  (see  Figure  3)  and  function 
Nje(x,y)  is  given  by  Equations  21  through  24. 

When  (x,  y)  is  on  a  boundary  segment  with  nodes  i  and  j  at  the  ends,  Equation  82 
is  simplified  to: 


fj(x,y)  =  fj,Nf  +f|jNJp  (83) 

where  Nf  =  (Sj  -  s)  /  Lp,  Nj  =  (s,  -  s)  /  Lp ,  and  s  is  the  relative  coordinate  along  P  and 
has  values  of  Sj  and  Sj  at  the  two  endpoints.  The  length  of  the  segment  is  LP  =  Isi  -  sjl. 

The  following  relations  referring  to  linear  function  Nie(x,y)  and  NiP(x,y)  are 
developed  for  later  use: 


9Ne  A 

VNf  =  M  + 
9  x 


9N*  - 

ay  J  " 


i  ( 


2Ae  V 


bi  i  +  Ci  j 


(84) 


r  r  Ae 

J  J  N'dxdy  =  y 


6Ae  /  60 

for  i  =  j  =  k 

|(  NfNjN^dxdy  =  < 

2Ae  /  60 

for  i  =  j  or  j  =  k  or  k  =  i 

e 

Ae  /  60 

for  i  j  *  k 

9NP  _  1 

9NP_  l 

9  s  Lp  ’ 

9  s  Lp 

(85) 


(86) 


(87) 


(88) 
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(89) 


J(Nr)JN'ds  =  jNf(N')Jds  =  lL' 

D  D 


(90) 


Now  substituting  the  Equations  80  and  81  into  Equation  78,  we  have 


+  -  JJbrjNfdaj  =  0  (i  =  1,...,N) 

(I)  (ID  (in) 


(91) 


where  e;  refers  to  those  elements  around  node  i  where  Nje(x,y)=0,  and  Pi  refers  to  the  two 
boundary  segments  on  each  side  of  node  i  when  i  is  a  boundary  node. 

For  element  e  with  three  nodes  i,  j,  and  k,  where  fj ,  given  by  Equation  82,  is 
substituted  into  the  i*,  j*  and  k*  Equation  of  91,  and  hence  the  second  terms  come  out 

n?  =  +  (J^aVNfVN'df^  +  (JJeaVNfVN*d£2ta 

=  Ayfj;  +  AyHj  +  Afkf)k  (92) 

nj  =  A'^+A^  +  A^  (93) 

nek  =  A^+A^  +  A^  (94) 


Likewise,  the  third  terms  in  Equation  91  become 
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mf  = 


-(JjebNf2do)fii  -({[BNrNJdO^  -(JJebNrNekdn)f\k 

=  B^f);  +  Byfjj  +  Bfkf)k 


(95) 


mj  =  Bft+BJV  BIkrik 


(96) 


mek=  B^+B^+B^ 


(97) 


Assuming  that  the  coefficients  a  and  b  also  vary  linearly  on  element  e,  i.e. 


a  =  SjN®  +  Sj  N®  +  ak  Nk 


b  =  bi  N*  +  bj  Nj  +  bk  Nk 


(99) 

(100) 


then  using  the  relations  (84),  (85)  and  (86),  we  have 


A'  =  JJi  VN'VNJdQ 

1  -(bjbj  +c,cJ)JJe|ai  N*  +  aj  N-  +  ak  Nek  jdn 

(I,J  =  i,j,k) 


4A 


e2 


If  'i 

=  ■^^■1  ai+  ai+  ak  +cicj) 


(101) 


Bu  =  -  JJJ  bi  Nf  +  bj  N-  +  bk  Nk 
A' 


N®NjdD 
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bi  +  bj+  bk+  2bi 


for  I  =  J 


Ae  f  ~  ~ 

bj  +  bj+  bk+  bi+  bj 


(I,  J  =  i,  j,  k) 


60 


for  I  *  J 


(102) 
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By  combining  the  expressions  for  r\  at  nodes  i,  j  and  k,  the  second  and  third  terms 
of  Equation  91  for  element  e  can  be  represented  by  means  of  matrix, 

[ne  +  nr]= ([ac]+ [b6])  {f|e}  ( e  =  1, 2, e  )  (103) 


where 


{f*C}  =  ^1* 


(104) 


and  the  expression  for  [Ae]  +  [Be  j  is  the  same  as  that  for  [K?  J ,  see  Equation  35. 


An  expression  for  the  first  term  of  Equation  91  may  be  obtained  by  applying 
certain  boundary  condition  and  also  using  Equations  87  through  90.  This  gives  the 
following  relationship  when  node  i  is  at  coastal  boundary. 


lCi 


=  -a J Nf^Si  Nf  +  aj  N[  j(Nf fjf  +  Nj’fjj’ )ds 


CP /i P  , 

ii^li  +QjTlj 


P  -  P 


(105) 


Similarly  for  node  j, 


icj = cjifir + cj^ 


p/^p 


(106) 


Therefore,  for  segment  P  with  node  i  and  j,  we  have  the  matrix  formula 


[icHc'FJ 


(107) 


where 
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(108) 


c  = 

'-'0  — 


-aj^ai  Nf  +  aj  Nj  JNfNfds 


'  Jj 

- afa  +  a 


(^  +  aj+2aI)  forI  =  J 


Lp  /„  „  x 

_T2  a  &1  +  Sj 


(I,  J  =  i,  j) 


for  I  *  J 


(109) 


When  node  i  is  at  open  boundary,  boundary  condition  (Equation  74)  applies,  and 
the  first  term  in  Equation  91  becomes 


In  =  J  a[^l  +  9  |qT  -  g]N.P  ds 


(HO) 


Similar  to  IcP  (Equations  107  and  108),  the  first  term  of  the  i*  and  j111  Equation  in 


(110)  is 


fc]-W} 


(111) 


where 


r,p  = 

*TT 


— p  ai+aj+2a,  forI  =  J 
12  v  / 


B  Tp 


l  r  _ 

TTP  %  +  aj 
12  v 


(I,  J  =  i,  j) 


for  I  *  J 


(112) 


The  second  terms  in  (1 10)  are 
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4 


J 


~  2 

1  N,ds  =  aqr 


as 


(n'M 

■as 

i  Jp  ds  ds  j 

and 


(113) 


~  2^f| 

aqr 


as 


(114) 


where  r  is  the  radius  of  the  semicircle.  These  terms  can  also  be  written  in  matrix  form  as 


KHr'FHD'] 


where 


IVp1_  aqr2  \l  -1 
l  2J  Lp  -1  1 


(115) 


(116) 


and 


mit 

asj 


(117) 


The  third  terms  in  Equation  1 10  are 

?„  =  -J  a  gNfds  =  -  a  gj  Nfds  =  -  j  a  gLp  =  if 

P  P  ^ 


therefore. 
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(118) 


[I?,Hr']=-?SgL' 


f\ 
1 


The  function  f\  will  be  obtained  by  first  computing  the  element  matrix  [Ae]and 

[Be]  for  elements  e  =  and  the  boundary  matrix  [cp],  [r,p],  [T2P],  [r3p]  for 

segments  P  =  1,...,  Np.  These  matrices  are  assembled  to  obtain  an  NxN  system  of 
equations, 


I([A']+[B-]){^}+X[cp}{^}+I([r;]+[r']}{f);}+I[D']+X[rn=0 

n  c  r  r  r 


Notice  that 


XK]= 


~ _ 2 


aqr 


9fi 
< - 

_df\ 

r  - . 

>  =  a  aqr 

'1  o' 

as 

A,’ 

AJ 

0  1_ 

{tv'HaJ7  (119) 


where  Ai  and  A2  are  two  points  that  connect  the  open  boundary  and  coastal  boundary  and 

the  wall  boundary  condition,  (9),  applies  to  df|/3s.  Therefore,  term  ^[Dp]  can  be 

r 

included  into  and  then  the  assembled  equation  becomes 

c 

Mf)MKJ]{f\cMK3]{lir}  =  {f}  020) 


or 

[A]{fi}={f} 


(121) 


This  linear  system  of  equations  may  be  solved  to  obtain  Tj. 
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4.4  Kinematic  Parameters 

Once  a  solution  for  f|  is  obtained  in  CGWAVE,  the  maximum  wave  velocity,  maximum 
wave  pressures,  wave  phase  angle  and  wave  amplitude  may  be  obtained  from  the  values 
of  f| .  These  quantities  are  obtained  as  follows. 

The  velocity  potential  for  water  particles  of  surface  water  waves  may  be  written  as 

<D(x,  y,  z,  t)  =  [<J>,  cos  cot  <j>2  sin  G)t]Z(z)  ( 1 22) 


where 


COSh[k(z  +  h)] 
cosh(kh) 


The  potential  may  be  expressed  in  terms  of  f|  by  substituting 


(123) 


,  g  „ 

O  =  —  T| 

ico 

into  Equation  122;  this  gives 

Z(z) 


®(x,  y,  z,t)  =— 9l| 
0)  1 


fl(x,y)e 


.(n  } 

-1  — KOt 

U  ) 


(124) 


(125) 


where 


=  (126) 

This  expression  may  be  written  as  follows  by  separating  the  real,  f)  j,  and  imaginary,  f)  2 
parts  on  f\ ,  and  replacing  -7t/2-cot  with  a;  this  gives 
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(127) 


0  =  —  [ri,  cos  a  +  r|2  sin  a]Z 
(0 


An  expression  for  the  velocity  of  water  particles  is  obtained  by  evaluating  the 
gradient  of  the  expression  for  <t>  in  the  last  equation;  this  gives 


v  =  *- 
x  CO 


(  ^7i 


LA 


dx 


cosa  + 


f  dr\:  ^ 

K  j 


sma 


g_ 

co 


dy 


IA 


cos  a  + 


fdrh) 

dy 


sin  or 


(128a) 


(128b) 


These  expressions  contain  the  horizontal  components  of  the  velocity.  For  simplicity,  Z  is 
taken  as  a  local  constant. 


The  magnitude  of  the  horizontal  components  of  the  velocity  is  obtained  by 
substituting  vx  and  vy  from  the  last  expression  into 

|v|2  =  (vx)2  +  (vy)2  (129) 


to  obtain 


lvf  = 


1  z2 

(G) 


Pn, 

dx 


+ 


l  dy 


V 


cos  a  + 


dn2 

dx 


2  f  Y 


+ 


dny 

dy 


sin2  a 


_  d%  df|2  [  dfj,  df|2 
dx  dx  dy  dy 

The  maximum  horizontal  velocity  occurs  at  locations  where  the  derivative  of  Ivl 
with  respect  to  a  is  equal  to  zero;  this  occurs  where 


sin(2a) 


(130) 


dfi,  XJ^L 


l  dx  J 


+ 


l  dy 


+ 


dfj2 

dx 


2  V 


dn, 

l  dy 


y  sin(2a) 
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Rearranging  terms  gives 


-2 


dr),  af|2  |  3fj,  df\2 

dx  dx  dy  dy 


cos(2a)  =  0 


(131) 


a  =  —  arctan 
2 


dfji  ajh  ,  ^2 

dx  dx  dy  dy 


ran,] 

2" 

(SSlY* 

f  dn2^ 

2“ 

V  dx  J 

dy  J 

l  dx  J 

dy  J 

(132) 


The  magnitude  of  the  horizontal  component  of  the  velocity  will  have  its  minimum  and 
maximum  values  at  this  value  of  a  and  at  a+n/2.  The  value  of  Ivl  is  calculated  at  both  of 
these  angles;  the  larger  value  is  the  maximum  velocity  over  all  times. 

The  pressure  is  obtained  from  the  linear  form  of  the  Bernoulli  equation; 

90  P 

- 1 - 1-  gz  =  constant  (133) 

at  p 

Note  that  the  pressures  associated  with  velocities  (dynamic  head,  V2V2)  are  ignored  in  this 
linear  form.  The  expression  for  <I>  in  Equation  125  is  substituted  into  this  expression  and 
terms  are  rearranged  to  obtain 

P  =  -pgz  +  pg9?(fje~1<0,)Z  +  constant  (134) 

The  maximum  pressure  over  a  wave  cycle  occurs  when  the  term  SR(f)e“,C0< )  is  equal  to 
H/2.  The  constant  is  chosen  such  that  the  hydrostatic  pressure  is  equal  to  zero  at  z=0; 
thus 


P™,=-Pgz  +  Pg"z  (135) 

The  wave  phase  angle  (5  is  obtained  from 
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(5  =  arctan  — 
Til 


(136) 


The  cosine  of  (3  varies  from  -1  to  1,  and  is  written  as  output. 

The  wave  amplitude  A  is  obtained  from 

A  =  -J%+  4  (137) 

A  ‘snap-shot’  of  the  sea  surface  elevation  at  time  =  0.0  is  obtained  from 

Tj  =  9t[f|e'i‘ot  j  =  [ti,  cos  cot  +  rj2  sin  cot]  (138) 


4.5  Floating  Docks 

A  floating  dock  inside  the  computational  domain  may  be  treated  following  the 
formulation  of  Tsay  &  Liu  (1983).  If  the  distance  from  the  bottom  of  the  bed  to  the 
surface  of  the  floating  body  is  d(x,y)  and  the  wave  number  corresponding  to  this  depth 
is  k(x,y) ,  the  governing  equation  for  region  under  the  floating  body  is 


V  •  (pVf))  =  0 


(139) 


where 


2 

fi  2kd  1 

1  + - = 

.k  j 

k  sinh  2kd , 

c2  =  gktanhkd 


(140) 

(141) 


Note  that  the  second  term  of  Equation  1  has  been  neglected  in  Equation  139.  This 
simplification  is  known  as  the  ‘rigid  lid’  approximation.  Tsay  and  Liu  (1983)  have 
demonstrated  that  this  formulation  gives  very  good  results  when  compared  with  some 
exact  solutions.  Equation  139  is  solved  in  CGWAVE  in  place  of  Equation  1  for  the 
region  where  a  floating  structure  is  present  in  the  domain  Cl. 
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5  ITERATIVE  CG  SOLUTION 


METHOD 


CGWAVE  provides  the  choice  of  choosing  the  wavelength  dependent  resolution 
to  solve  the  governing  equation.  An  in-depth  investigation  was  performed  to  determine 
the  sensitivity  of  the  solution  as  a  function  of  wavelength.  Based  on  our  analysis,  it  is 
recommended  that  ten  or  more  points  per  wavelength  to  be  used.  In  regions  where  the 
change  in  topography  or  wave-amplitude  is  rapid,  a  finer  grid  is  necessary.  To  solve  the 
elliptic  mild-slope  equation  without  approximating  the  physics  (as  in  parabolic  models), 
the  problem  has  to  be  solved  simultaneously  over  the  entire  domain.  This  leads  to  a  very 
large  linear  system  of  equations  (e.g.  Equation  69).  Direct  methods  (e.g.  Gaussian 
elimination)  are  often  inapplicable  due  to  extremely  large  storage  requirement  of  the 
matrix  [A],  Iterative  methods,  on  the  other  hand,  require  memory  for  only  the  non-zero 
elements  in  [A].  Since  [A]  is  highly  sparse,  iterative  methods  can  significantly  enhance 
the  ability  of  the  elliptic  type  mild-slope  equation  models  to  handle  large  domain 
problems.  However,  most  iterative  procedures  require  [A]  to  be  diagonally-dominant  or 
symmetric  and  positive-definite.  Unfortunately,  the  coefficient  matrix  [A]  is  not 
diagonally-dominant,  nor  symmetric  and  positive  definite.  Panchang  et  al.  (1991) 
recommended  that  the  following  Gauss  transformation  is  applied  to  Equation  69  : 

[A'jA]{fi}=[A']{f}  (142) 


where  A*  is  the  complex  conjugate  transpose  of  [A].  Xu  (1995)  has  shown  that  the  new 
coefficient  matrix  [A*]  [A]  is  Hermitian  and  positive-definite.  Therefore,  the  modified 
conjugate-gradient  method,  which  is  often  several  orders  of  magnitude  faster  than  many 
other  schemes,  including  the  traditional  conjugate  gradient  scheme,  is  guaranteed  to 
converge  when  applied  to  Equation  125.  The  algorithm  is  implemented  in  CGWAVE  as 
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follows: 


1.  Select  trial  values  for  { fj  0}  (i.e.  0th  iteration)  for  all  nodes  in  model  domain 
where  the  solution  is  desired. 

2.  Compute  for  all  points  the  residual  {r0}  =  {f}  -  [A]  { fj  0}and  the  left  hand  side 
of  Equation  125  as  {p0}  =  [A*]{ro}. 

th  ||[A*]{ri}||2 

3.  Compute  for  the  i  iteration  the  parameters  a*,  defined  as:  <X:  = 11 - V 

|[A](P,)f 

4.  Update  {Tii+1}  =  {fjj  +  ajlpj}  for  all  points. 

5.  Check  for  convergence  of  the  solution.  The  criterion  for  convergence  used  in 
CGWAVE  is 


|[A]|f|M}-|f}f 

IKM! 


(143) 


where  8  is  a  prescribed  tolerance.  If  Equation  126  is  satisfied,  stop. 

6.  If  Equation  126  is  not  satisfied,  compute,  for  each  grid  point, 
{ri+1 }  =  {v,}  -  ccJAHPi}. 


7.  Compute  for  the  i*  iteration 


P, 


|A-j{rM)f 
||[A-]{r,}f  ' 


8.  Compute  {pi+, }  =  [A"  ]  {rj+I )  +  P,  (p, } . 

9.  Set  i  =  i  +  1,  and  go  to  step  3. 


In  the  procedure  above,  the  module  of  an  array  ||{x}||  is  defined  as 


This  procedure  has  been  demonstrated  to  work  very  well  in  finite-element  models. 
In  fact,  these  techniques  were  found  to  be  far  more  efficient  in  the  present  finite-element 
runs  than  for  the  finite-difference  studies  of  Xu  &  Panchang  (1993)  and  Panchang  et  al. 
(1991). 


Li  (1994)  has  recently  suggested  modifications  to  the  iterative  schemes  of 
Panchang  et  al.  (1991)  for  enhancing  the  convergence.  Comparison  between  these  two 
schemes  were  given  in  Xu  and  Panchang  (1995).  Although  these  modifications  resulted  in 
faster  convergence,  it  was  noticed  that,  unlike  the  basic  procedures  of  Panchang  et  al. 
(1991),  they  do  not  converge  to  a  solution  monotonically.  Instead,  the  residual  error 
decreases  in  an  oscillatory  manner  as  the  iterations  proceed.  Also,  the  gain  in  CPU  time 
varies  from  case  to  case  and  appears  to  be  problem-specific,  but  it  can  save  over  50%  of 
CPU  time  for  some  applications.  This  scheme  has  also  been  included  in  CGWAVE  as  an 
alternative  choice,  and  consequently,  there  are  two  types  of  solvers  provided  in 
CGWAVE. 

The  solving  algorithm  is  slightly  different  when  non-linear  mechanisms  are 
incorporated  into  the  solution.  In  this  case,  the  solution  is  first  solved  as  if  there  were  no 
non-linear  mechanisms.  This  solution  is  then  used  to  prescribe  initial  conditions  for  the 
non-linear  mechanisms,  and  the  system  of  equations  are  modified  to  incorporate  the 
resulting  non-linear  mechanisms.  CGWAVE  then  solves  using  the  modified  equations. 
This  iterative  method  of  solving,  modifying  the  system  of  equations  with  non-linear 
mechanisms,  and  solving  again  continues  until  the  resulting  solution  no  longer  changes 
between  iterations. 

The  solving  algorithm  is  also  different  when  spectral  wave  conditions  are  used.  In 
this  case,  a  solution  is  first  obtained  for  each  individual  spectral  component.  The  final 
solution  is  obtained  by  linear  superposition  of  the  solutions  for  all  spectral  components. 
This  superposition  is  performed  after  a  spectral  run  is  completed,  using  a  post  processing 
program. 
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6  GENERATION  OF 


FINITE-ELEMENT  NETWORK 


CGWAVE  requires  a  two-dimensional  (2D)  triangular  grid  network  for  its  finite- 
element  calculations.  Although  several  grid  generation  packages  are  available,  they  are 
not  suitable  for  elliptic  coastal  wave  models  for  which  the  size  of  the  elements  must  be 
related  to  the  wave  length  (which  varies  with  local  water  depth)  for  proper  resolution.  A 
semicircular  open  boundary  has  to  be  created  for  special  open  boundary  treatment,  and 
reflection  coefficients,  which  may  vary  from  one  part  of  the  coastal  boundary  to  another, 
are  also  required  as  input  data  for  the  model.  To  deal  with  these  special  problems, 
CGWAVE  has  been  interfaced  with  the  grid-generator  associated  with  the  SMS  (Surface 
water  Modeling  Systems)  flow  modeling  package.  The  Engineering  Computer  Graphics 
Laboratory  is  developing  this  state-of-the-art  package  for  the  US  Army  Corps  of 
Engineers  at  Brigham  Young  University. 

SMS  contains  a  set  of  2D  hydrodynamic  models  and  a  general  purpose  grid 
generation  and  visualization  package.  SMS  includes  an  efficient  finite-element  grid- 
generator.  However,  this  grid-generator  was  originally  designed  for  other  types  of 
hydrodynamic  models.  Three  utility  programs  that  help  interface  CGWAVE  with  the 
SMS  grid-generator  have  been  developed  for  use  outside  SMS,  prior  to  the  full 
integration  of  CGWAVE  into  SMS.  Given  a  coarse  rectangular  array  of  bathymetric  data, 
these  programs  generate  a  wavelength-dependent  triangular  nodal  network  (based  on  the 
user-specified  resolution,  i.e.  the  number  of  points  per  wave  length),  automatically 
construct  the  semi-circular  open  boundary,  assign  reflection  coefficients  along  the  coastal 
boundaries,  eliminate  unwanted  land  points,  etc.  The  resulting  grid  and  boundary  data 
from  SMS  are  then  filtered  by  another  utility  program  for  use  by  the  wave  model.  The 
output  from  the  wave  model  can  be  processed  and  then  plotted  by  using  SMS.  This  makes 
model  implementation  very  efficient  and  allows  the  user  to  view  a  graphic  representation 
of  the  solution. 
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7  PROCEDURE  FOR  STAND-ALONE 


USAGE  OF  CGWAVE 


CGWAVE  is  an  efficient  and  easy-to-use  wave  model.  Although  more  advanced 
features  are  anticipated  to  evolve  in  the  future,  the  current  version  is  sufficiently  user- 
friendly  and  can  provide  the  base  information  for  shallow  water  wave  propagation  and 
transformation  for  harbor/coastal  problems  or  offshore  open  sea  problems.  The 
information  provided  by  this  model  can  be  further  used  to  calculate  wave  forces  on 
structures,  stability  of  breakwaters,  sediment/pollutant  transport  and  can  be  used  to  assist 
harbor/coastal  protection  structure  design. 

In  this  section  of  the  manual,  step-by-step  instructions  are  given  for  use  of 
CGWAVE.  The  procedure  outlined  here  pertains  to  the  usage  of  CGWAVE  as  stand¬ 
alone  (i.e.  outside  SMS);  this  usage  is  referred  to  as  MODE-1.  MODE-1  usage  does 
require  the  use  of  SMS  for  some  procedures.  A  similar,  automated  procedure  is  available 
within  SMS,  and  is  referred  to  as  MODE-2.  After  successful  installation  of  CGWAVE 
and  SMS,  the  user  can  proceed  with  applications  of  CGWAVE.  In  general,  the  MODE-1 
application  of  CGWAVE  consists  of  eight  steps: 

Step  1.  After  a  region  for  simulation  has  been  chosen,  a  file  containing  (x,y,z 
where  z  =  depth)  data  in  rectangular  fashion  may  be  created  from  a  map  or  other  data 
base.  The  input  file  format  is  described  in  detail  in  section  7.2.2  of  this  manual  and  in  the 
File  Formats  section  of  the  SMS  Reference  Manual. 

Step  2.  The  raw  data  gathered  in  Step  1  is  processed  by  the  program  “resol”. 
Coastal  and  open  boundaries  are  defined  and  most  extraneous  data  is  eliminated. 
Wavelength  dependent  resolution  is  established,  and  the  revised  (x,y,z)  data  set  is 
presented  for  constructing  the  finite  element  grid. 

Step  3.  The  finite  element  grid  is  constructed  using  SMS,  and  boundary 
conditions  are  assigned  to  the  open  and  coastal  boundaries. 
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Step.  4.  CGWAVE  input  parameters  are  specified.  These  parameters  include: 
incident  wave  condition,  nonlinear  mechanisms  (wave  breaking,  bottom  friction  and 
nonlinear  dispersion),  choice  of  solvers,  choice  of  open  boundaiy  condition,  iteration 
control  parameters,  etc. 

Step  5.  All  relative  information  needed  for  input  to  the  model  is  consolidated  into 
a  single  input  file.  This  processing  is  done  by  program  “reform”. 

Step  6.  Run  the  wave  model  and  obtain  an  output  wave  potential  file. 

Step  7.  Process  the  output  file  (using  program  “trans”)  to  obtain  amplitude,  phase 
and  kinematics  solution  files. 

Step  8.  Use  SMS  to  view  various  solution  images  using  a  variety  of  contouring 
options. 

7.1  File  Descriptions 

The  following  files  contain  source  code  for  the  CGWAVE  modeling  package: 

cgwave.f  —  Fortran  source  code  of  CGWAVE. 
cgwave.gen  —  Template  of  the  general  information  file, 
reform.f  —  Fortran  source  code  of  reform, 

resol.f  —  Fortran  source  code  of  resol. 

trans.f  —  Fortran  source  code  of  trans. 

Given  below  is  a  list  of  example  files  for  application  of  CGWAVE.  Although  one 
can  use  any  name  for  them,  the  following  convention  is  used  throughout  this  manual.  It  is 
also  recommended  that  the  user  follow  this  convention  for  simplicity  and  consistency. 

filename.xyz  —  File  to  be  imported  to  SMS  and  to  be  used  by  program  “resol”. 

Contains  (x,y, depth)  data  describing  the  bathymetry. 
filename.be  —  Boundary  condition  file  created  by  SMS. 
filename.geo  —  Geometry  file  created  by  SMS. 
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filename.gen  —  General  information  file,  text  format.  Contains  required  input 
control  parameters  forCGWAVE. 

filename.dat  -  CGWAVE  and  “trans”  input  data  file.  Created  by  program 
“reform”. 

filename.out  -  CGWAVE  wave  potential  output  file.  Input  to  program  “trans” 
filename.txt  ~  Text  output  file  created  by  trans.  Contains  general  information 
about  the  project  modeled,  CGWAVE  input  parameters,  solution 
status  and  complete  output  of  x,y  depth,  amplitude  and  phase  at 
each  node  in  text  format. 

filename.sol  -  Solution  output  file  created  by  program  “trans”.  Contains 
information  about  computed  amplitude,  phase  and  kinematics. 
May  be  read  by  SMS  for  graphic  presentation  of  the  solution. 
ASCH  text  format. 

filename.pls  -  User  created  input  file  for  program  “trans”  which  provides 
information  about  selected  node  ids,  points,  lines 
and  rectangular  sub-regions  where  output  is  desired, 
filename.sel  —  Selected  locations  output  file  created  by  program  “trans”. 
filename.flw  -  Floating  dock  solution  output  created  by  program  “trans”.  Of 
relevance  only  if  there  is  a  floating  dock  inside  the  computational 
domain.  Velocity  potential  output  is  available  in  this  sub-domain, 
which  may  be  further  used  (outside  CGWAVE)  to  calculate  forces 
and  moments. 

filename.alp  —  Animated  film  loop  for  view  of  variations  in  amplitude  and  phase 
with  time.  For  use  in  SMS. 

7.2  Raw  Data  Processing  -  Program  “resol” 

Grid  resolution  is  dependent  on  the  local  water  depth.  In  general,  10  or  more 
nodes  per  (local)  wavelength  are  desired.  In  most  cases,  a  large  number  of  extra  nodal 
points  must  be  created  over  the  raw  data  points  to  satisfy  the  resolution  requirements. 
Hence  it  is  desirable  to  have  these  points  generated  before  using  SMS  to  create  a  finite 
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element  network.  Nodes  with  negative  or  zero  depths  (dry  land)  are  prohibited  within  the 
computational  domain.  A  minimum  water  depth  may  be  judicially  chosen  and  regions 
under  this  minimum  depth  (higher  in  elevation)  should  be  eliminated.  For  the  open 
boundary,  it  is  necessary  to  have  properly  distributed  points  created  on  a  semicircle  (or 
full  circle  for  offshore  applications),  and  points  outside  the  open  boundary  must  be 
discarded.  These  tasks  are  theoretically  trivial,  but  very  time  consuming.  The  utility 
program  “resol”  accomplishes  these  tasks  quickly  and  with  minimal  effort  on  the  part  of 
the  user. 


7.2.1  Viewing  XYZ  Data  and  Selecting  Open 
Boundary  End  Points 

Consider  as  an  example  the  raw  data-set  called  “swansr.xyz”.  This  raw  data  may 
be  viewed  in  SMS.  This  can  be  done  as  follows.  Run  the  SMS  executable  file.  When 
SMS  is  open,  a  tool  palette,  graphics  window,  and  pull-down  menus  at  the  top  of  the 
screen  are  visible.  From  the  “File”  menu  choose  “import”.  Select  “XYZ  data”  from  the 
“Select  Import  Format”  window.  Select  the  drive,  directory  and  file  “swansr.xyz”  using 
the  “Open”  window.  From  the  tool  palette  select  the  “Create  Linear  Triangle”  tool  Note 
that  the  name  of  the  tool  that  is  currently  under  the  mouse  pointer  is  displayed  at  the 
bottom  of  the  SMS  window.  From  the  “Elements”  menu  select  “Triangulate”.  This 
results  in  a  finite  element  triangular  network.  Select  the  “Display  Options”  tool  from 
either  the  tool  palette  or  from  the  “Display”  pull-down  menu.  The  “Display  Options” 
window  allows  the  user  to  view  or  hide  various  display  features.  Toggle  off  “Nodes”. 
Toggle  the  “Contours”  on.  Click  on  the  “Options”  button  to  the  right  of  the  “Contours” 
toggle.  Select  “Color  fill  between  contours”,  click  “OK”.  Click  “OK”  again  to  exit  the 
“Display  Options”  window.  A  color  representation  of  the  triangulated  grid  is  now  visible. 
It  is  suggested  that  the  user  become  familiar  with  the  function  of  the  tools  and  options 
before  attempting  any  serious  modeling  effort.  Before  exiting  SMS,  the  user  should  write 
down  some  information  that  is  used  to  create  the  semicircular  open  boundary.  Pick  a 
minimum  water  depth  and  roughly  estimate  the  starting  and  end  point 
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coordinates  (counter-clockwise)  of  the  semicircular  open  boundary.  Since  the  wave 
model  does  not  allow  zero  or  negative  depths  in  the  computational  domain,  a  minimum 
depth  must  be  specified  when  executing  the  program  “resol”.  Make  sure  that  both 
semicircle  end  points  are  inside  an  area  where  water  depth  is  greater  than  the  minimum 
depth  on  at  least  one  node  of  the  background  “rectangle”.  In  this  example,  2.0  (meters)  is 
used  as  the  minimum  depth  and  the  points  (3500.0,2700.0),( 100.0,2450.0)  as  the  starting 
and  end  point  coordinates  for  the  semicircle.  You  may  now  quit  the  SMS  (do  not  save 
any  file  information  at  this  time). 

7.2.2  XYZ  File  Format 

The  input  data  file  for  the  program  “resol”  must  be  a  rectangular  domain  with 
constant  dx  and  dy  for  each  direction.  The  file  has  the  following  format: 

XYZ 
xl  yl  zl 
x2  y2  z2 


xn  yn  zn 


For  example,  by  viewing  “swansr.xyz”  with  a  text  editor,  the  following  will  be  seen: 


XYZ 


0.0 

0.0 

0.0 

0.0 

200.0 

0.0 

0.0 

400.0 

0.0 

0.0 

600.0 

0.0 

0.0 

800.0 

0.0 

0.0 

1000.0 

0.0 

0.0 

1200.0 

0.0 

0.0 

1400.0 

0.0 

0.0 

1600.0 

0.0 

0.0 

1800.0 

0.0 

0.0 

2000.0 

0.0 

0.0 

2200.0 

0.0 
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0.0 

2400.0 

0.0 

0.0 

2600.0 

15.3 

0.0 

2800.0 

15.3 

0.0 

3000.0 

15.3 

0.0 

3200.0 

14.9 

0.0 

3400.0 

23.8 

0.0 

3600.0 

23.8 

0.0 

3800.0 

27.5 

0.0 

4000.0 

15.3 

0.0 

4200.0 

24.4 

0.0 

4400.0 

13.7 

200.0 

0.0 

0.0 

200.0 

200.0 

0.0 

200.0 

400.0 

0.0 

Note  that  this  is  only  the  top  portion  of  the  swansr.xyz  file.  The  first  line  must  be  “XYZ” 
followed  by  (x  y  z)  data.  This  xyz  file  should  be  oriented  such  that  the  x-axis  points  to 
the  east  on  the  screen  when  the  data  set  is  displayed. 


7.2.3  Program  “resol” 

Execution  of  program  “resol”  will  result  in  the  creation  of  a  new  XYZ  data  file 
called  “swanst.xyz”,  which  is  to  be  used  as  input  to  SMS  for  further  grid  generation 
operations.  This  file  has  better  resolution  than  the  raw  data  file  (e.g.  swansr.xyz),  and 
contains  information  regarding  the  semicircular  open  boundary. 

To  begin  the  example,  run  the  “resol”  executable.  A  program  information  header 
will  be  displayed  on  the  screen.  The  following  information  will  appear  after  the  header: 

Input  filename  :  (*.xyz) 

swansr.xyz  (Hit  <Enter>  after  each  response) 

Output  Filename:  (*.xyz) 

swanst.xyz 

Input  data  spacing. 

Dx  is: 

200.0 

Dy  is: 
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200.0 


Incident  wave  period 
30.0 

Enter  type  of  open  boundary: 

Full  Circle  (=0),  Semi  Circle  (=1) 

1 

Start  from: 

3500.0  2700.0 

End  at  (ccw): 

100.0  2450.0 

Are  they  Fixed  (=0)  or  Approximate  (=1)  ? 

1 

Number  of  Nodes  per  Wavelength  : 

10 

Interpolate  between  land  and  water  ? 

(yes=0,  no=l) 

0 

Minimum  Water  Depth : 

2.0 

Minimum  Node  Spacing : 

.001 

Starting  point  of  the  semicircle  : 

3511.111  2711.111 

End  point  of  the  semicircle  : 

80.000  2440.000 

End  of  program.  Press  Enter  to  continue 
Some  notes  regarding  the  above  questions  and  answers  follow. 


Minimum  Water  Depth.  When  specifying  a  minimum  water  depth,  the  estimated 
semicircle  end  points  may  or  may  not  satisfy  the  minimum  depth  criteria.  In  most  cases. 
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these  coordinates  will  be  either  inside  or  outside  the  domain  (and  rarely  on  the  exact 
boundary  indicated  by  the  minimum  depth.)  By  choosing  “Approximate”,  the  program 
“resol”  will  locate  a  boundary  point  closest  to  the  point  specified.  Note  also  that  program 
“resol”  requires  that  the  location  of  the  estimated  points  should  lie  inside  a  rectangle  (of 
the  background  or  raw  data  network)  with  at  least  one  comer  having  a  depth  greater  than 
or  equal  to  the  minimum  depth. 

Number  of  Nodes  per  Wavelength.  It  is  strongly  recommended  that  at  least  10  or  more 
points  per  wavelength  be  used  whenever  possible.  A  resolution  of  7  or  8  points  per 
wavelength  may  yield  a  reasonable  solution,  but  6  should  be  treated  as  an  absolute 
minimum.  15  points  per  wavelength  are  sufficient  for  most  applications. 

Interpolate  Between  Land  and  Water.  Other  than  some  ideal  cases  with  constant  depth 
(e.g.  laboratory  or  analytical  test  cases),  always  answer  yes  to  this  question. 

Minimum  Node  Spacing.  This  feature  allows  the  user  to  limit  the  density  of  nodes  (and, 
indirectly,  the  total  number  of  nodes)  placed  within  the  domain.  This  may  be  necessary 
for  extremely  shallow  regions  or  for  very  large  domains.  This  feature  can  contradict  the 
“Number  of  Nodes  per  Wavelength”  parameter,  and  should  be  used  with  caution.  It  is 
recommended  that  the  user  set  this  parameter  to  a  very  small  value  (.001  or  smaller)  for 
most  applications. 

7.3  Creation  of  Finite  Element  Grid 

A  description  of  the  steps  involved  in  creating  a  finite  element  grid  and 
constructing  CGWAVE  input  files  using  SMS  follows.  This  procedure  describes 
MODE-1  usage  of  CGWAVE.  It  is  highly  recommended  that  the  user  be  fully  familiar 
with  the  MODE-1  usage,  since  the  automated  procedure  for  MODE-2  makes  use  of  the 
steps  described  below.  In  MODE-1,  the  user  needs  to  create  the  following  files: 

*.geo  (geometry  file) 
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*.bc  (boundary  condition  file) 

*.gen  (general  input  file) 

These  files  are  used  by  the  program  “reform”  to  construct  the  main  CGWAVE  input  data 
file.  Specific  steps  pertinent  to  creating  each  of  the  above  files  are  provided  below. 

7.3.1  Geometry 

Open  SMS.  From  the  “File”  menu  choose  “import”.  Select  “XYZ  data”  from  the 
“Select  Import  Format”  window.  Select  the  drive,  directory  and  file  “swanst.xyz”  in  the 
“Open”  window.  From  the  tool  palette  select  the  “Create  Linear  Triangle”  tool.  From  the 
“Elements”  menu  select  “Triangulate”.  Any  nodes  and  elements  that  lie  outside  the 
desired  model  domain  must  be  manually  removed.  Using  the  “Select  Element”  and 
“Select  Nodes”  tools  does  this.  The  <Shift>  and  <Ctrl>  keys  may  be  used  in  combination 
with  these  tools  to  remove  multiple  and  groups  of  elements  or  nodes.  For  example,  a 
group  of  elements  may  be  selected  for  deletion  as  follows:  Hold  the  <Ctrl>  key,  click  (do 
not  release)  the  mouse  button  on  the  center  of  an  element.  Drag  the  pointer  across  the 
group  of  elements  to  be  deleted.  An  arrow  beginning  at  the  first  selection  point  should 
follow  the  pointer.  Any  element  touched  by  this  arrow  will  be  selected.  Elements  may  be 
added  to  (or  removed  from)  a  selection  group  by  holding  the  <Shift>  key  while  selecting 
an  element.  The  selected  elements  are  then  deleted  with  the  <Del>  key.  Other  selection 
methods  are  available  in  the  “Edit”  and  “Elements”  menus.  Remember  to  clean  up  and/or 
remove  the  elements  in  areas  where  islands  or  structures  are  located.  After  obtaining  a 
satisfactory  grid-network,  boundary  conditions  must  be  specified,  as  detailed  in  the 
following  section.  Select  “Save  Geometry”  in  the  “RMA2”  menu.  Name  the  file 
swanst.geo. 


54 


7.3.2  Boundary  Conditions 

Defining  Open  and  Coastal  Boundaries.  The  boundary  conditions  for  this  grid 
may  now  be  specified.  Pick  “Create  Nodestring”  tool  from  the  tool  palette.  This  tool  is 
used  to  identify  the  nodes  along  the  edge  of  the  model  domain.  Use  the  mouse  pointer  to 
select  the  start  point  (node)  of  the  semicircular  open  boundary.  (The  “Zoom”  tool  may  be 
used  to  help  identify  the  start  node.)  Press  and  hold  the  <Ctrl>  key  and  select  the  end 
point  (node)  of  the  semicircular  open  boundary.  Press  <Enter>.  The  open  boundary  has 
been  defined  by  a  single  “nodestring”.  For  open  sea/offshore  problems,  where  open 
boundary  is  a  full  circle,  similar  procedure  applies.  For  this  example,  we  will  next  define 
the  coastline  boundary  using  a  single  string.  Proceed  with  the  above  selection  procedure, 
this  time  using  the  end  point  of  the  semicircle  as  the  new  start  point.  The  new  end  point 
is  the  start  point  of  the  semicircle.  In  some  modeling  situations,  the  user  may  want  to 
assign  different  reflection  coefficients  to  sections  of  the  coastline.  If  this  is  the  case,  a 
nodestring  must  be  created  for  each  area  of  differing  reflection  coefficient.  A  new 
nodestring  must  always  begin  with  the  end  point  of  the  previous  string.  The  direction  of 
nodestrings  should  be  specified  such  that  when  moving  from  the  start  point  to  the  end 
point  of  the  nodestring,  the  modeled  domain  is  on  the  left  side  of  the  nodestring,  i.e. 
counter-clockwise. 

Defining  Islands  Boundaries.  Island  boundaries  are  defined  in  much  the  same 
way  as  described  above,  with  the  following  exceptions.  Island  nodestrings  are  created 
such  that  the  modeled  domain  is  on  the  right  side  of  the  nodestring  when  moving  from  the 
start  point  to  the  end  point,  i.e.  counter-clockwise.  A  complete  island  boundary  must 
consist  of  two  or  more  nodestrings,  as  SMS  does  not  allow  a  single  nodestring  to  begin 
and  end  on  the  same  node.  Islands  should  be  defined  after  all  open  and  coastal 
boundaries  have  been  defined. 

Assigning  Reflection  Coefficients.  Pick  “Select  Nodestring”  tool.  A  small 
selection  area  (box)  should  appear  in  the  center  of  each  nodestring.  Place  the  pointer  on 
the  selection  box  for  the  open  boundary  and  select  it.  Go  to  the  “RMA2”  menu  and  select 
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“Assign  BC...”.  In  the  “RMA2  Assign  Boundary  Condition”  window,  pick  the  ‘Head 
BC”  option.  Set  the  elevation  (lower  part  of  the  window)  to  5.0.  This  value  is  used  to 
tell  CGWAVE  that  this  nodestring  is  an  open  boundary.  Click  “OK”.  Repeat  this 
procedure  for  all  other  boundaries.  For  all  boundaries  other  than  the  open  boundary,  the 
elevation  value  represents  the  reflection  coefficient  and  should  be  between  0.0  and  1.0, 
corresponding  to  the  reflection  coefficient  of  that  portion  of  coastal  boundary  (including 
islands).  The  value  of  1.0  for  reflection  coefficient  should  be  used  for  a  fully  reflecting 
boundary  and  0.0  used  for  fully  absorbing  boundary.  Save  the  boundary  condition 
information  by  selecting  “Save  BC...”  from  the  “RMA2”  menu.  For  the  example  case, 
name  the  file  “swanst.bc”. 

Floating  Structures.  Floating  structures  are  discussed  in  detail  in  section  4.5. 
When  a  floating  structure  is  included  inside  the  domain,  the  elements  containing  the 
floating  structure  must  be  both  marked  and  adjusted  to  signify  and  account  for  the 
presence  of  the  structure.  The  elements  contained  by  the  area  of  the  floating  structure  are 
assigned  a  “Material  Type  ID”  of  “3”  as  follows.  First,  select  all  the  elements  within  an 
area  bounded  by  the  perimeter  of  the  floating  structure.  Then  select  “Assign  Material 
Type”  from  the  “Elements”  menu.  Click  on  the  “New”  button  in  the  “Materials  Editor” 
until  the  “ID:”  changes  to  “3”.  Click  “Select”.  Next,  modify  the  depth  of  the  nodes  in 
this  area  as  described  in  section  4.5.  The  new  depths  are  shown  in  Figure  6  as  d .  This 
can  be  done  using  the  “Transform  Mesh”  option  found  in  the  “Nodes”  menu.  CGWAVE 
will  now  recognize  this  area  as  a  floating  structure  and  treat  it  accordingly. 

Renumbering  the  Geometry.  After  completing  the  grid  generation  and  boundary 
definition,  the  grid  must  be  renumbered.  Using  the  “Select  nodestring”  tool,  select  any  of 
the  nodestrings.  From  the  “Elements”  menu  choose  “Renumber”.  In  the  “Renumbering 
Opts”  window,  select  “Band  Width”  and  then  “OK”.  After  renumbering,  the  grid 
generation  process  is  complete.  Re-save  the  geometry  and  boundary  condition  files  using 
the  same  names  given  before  (swanst.geo  and  swanst.bc).  Some  basic  information  about 


the  grid  may  be  viewed  by  clicking  the  “Display  Module  Info”  tool,  or  by  selecting  “Get 
Info...”  from  the  “File”  menu.  You  may  now  quit  SMS. 

7.4  The  General  Information  File  (*.gen) 

A  template  called  “cgwave.gen”  is  provided  for  the  benefit  of  the  user.  A  new 
copy  of  this  file  should  be  made  for  each  application.  The  contents  of  the  file  are  written 
in  ASCH  text  and  appear  below.  Note  that  the  incident  wave  conditions  contain  data  for  a 
single  (monochromatic)  model  run.  Spectral  input  conditions  may  be  input  by  including 
additional  lines  below  the  first  (Angle/Period/Amplitude)  line.  Two  empty  (blank)  lines 
should  follow  the  last  (Angle/Period/Amplitude)  line,  separating  it  from  the  “2)  Open 
Boundary  ...”  line.  Incident  angles  are  measured  in  degrees  and  the  incident  direction  is 
defined  such  that  the  incident  wave  travels  in  the  positive  x-direction  when  0;  is  equal  to 
zero. 

PROJECT :  Project  Name 
(...blank  line) 

(...blank  line) 

1)  Incident  Wave  Conditions: 

(...blank  line) 

No.  Incident  Wave  Angle  II  Wave  Period  (s)  II  Incident  Wave  Amplitude 
1  0.0  30.0  1.0 

(...blank  line) 

(...blank  line) 

2)  Open  Boundary  (Bessel  Series  =  0,  Parabolic  =  1,  Box  =  2) :  1 
(...blank  line) 

If  Bessel  Series,  Number  of  Terms  Used  :  50 

If  Parabolic,  Exterior  Coastline  Reflection  :  0.0 
(...blank  line) 

(...blank  line) 

3)  Bottom  Friction  (Yes=l,No=0)  ?  0 
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If  Yes,  Friction  Coefficient :  0.12 
(...blank  line) 

(...blank  line) 

4)  Wave  Breaking  (Yes=l,No=0)  ?0 
(...blank  line) 

If  Yes,  Breaking  Coefficient :  0. 15 
(...blank  line) 

(...blank  line) 

5)  Nonlinear  Dispersion  Relation  (Yes=l,No=0)  ?  0 
(...blank  line) 

(...blank  line) 

6)  Solver  (eg  standard  =  0,  modification  =  1) :  1 
(...blank  line) 

Tolerance  of  Convergence  for  Linear  Equations  :  0.1  E-8 
Interval  for  Checking  Convergence  :  100 
Maximum  Iterations  for  the  Linear  Equations  :  100000 
(...blank  line) 

(...blank  line) 

7)  Tolerance  of  Convergence  for  Nonlinear  Mechanisms  :  0.0001 
(...blank  line) 

Maximum  Iterations  for  the  Nonlinear  Mechanisms  :  8 
(...blank  line) 

(...blank  line) 
end  of  data 

This  file  provides  the  basic  incident  wave  and  boundary  type  information  and 
various  control  parameters  for  CGWAVE.  Due  to  the  programming  methods  used  to  read 
the  data  from  this  file,  the  text  must  not  be  altered.  The  parameter  values  may  be 


replaced  with  values  that  are  specific  for  the  new  application.  All  the  numbers  must  be 
present  and  in  the  correct  location,  whether  they  are  used  or  not. 

7.5  Assembly  of  Input  Information  -  Program  “reform” 

The  program  “reform”  is  used  for  assembling  and  converting  all  input  information 
into  the  format  needed  by  CGWAVE.  It  reads  the  geometry,  boundary  condition  and 
general  information  files  and  constructs  a  unified  input  data  file  for  CGWAVE.  The 
program  also  verifies  some  components  of  the  provided  data.  Program  “reform”  makes 
the  central  part  of  CGWAVE  portable.  In  other  words,  a  workstation  may  be  used  to 
create  and  manipulate  the  graphics,  and  CGWAVE  may  be  run  on  any  other  platform 
(PC,  mainframe,  etc.).  Also,  because  “reform”  sorts  and  rewrites  the  data  in  clean,  nested 
groups,  the  same  data  set  can  be  used  to  perform  simulations  with  different  incident  wave 
directions,  different  combinations  of  boundary  reflection  coefficients,  different  incident 
wave  amplitudes  and  periods  (to  some  extent,  since  resolution  is  wave  period  dependent), 
etc.  without  going  through  the  trouble  of  editing  the  general  information  file  and  re¬ 
running  “reform”  for  each  case.  This  is  done  by  simply  editing  a  few  lines  in  the  top 
sections  of  the  data  (*.dat)  file.  Another  feature  of  “reform”  that  is  especially  useful  is 
the  ability  to  check  the  grid  resolution  for  the  overall  domain. 

Here,  the  rectangular  harbor  test  is  used  to  demonstrate  of  the  function  of 

“reform”.  Run  the  program  “reform”.  After  the  informational  header,  the  following 

appears  on  the  screen: 

General  information  filename:  (*.gen) 
rect.gen 

Boundary  condition  filename  from  SMS:  (*.bc) 
rect.bc 

Geometry  filename  from  SMS:  (*.geo) 
rect.geo 

Output  data  filename  for  CGWAVE:  (*.dat) 
rect.dat 
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Do  you  want  resolution  analysis  (yes/no)  ? 
yes 

Number  of  points  per  wave  length  distribution  (%) 


From  0  to  5 

.0 

From  5  to  6 

.0 

From  6  to  7 

.0 

From  7  to  8 

.0 

From  8  to  9 

.0 

From  9  to  10 

.0 

From  10  to  11 

.0 

From  11  to  12 

.0 

From  12  to  13 

.0 

From  13  to  14 

.0 

From  14  to  15 

.0 

From  15  up 

100.0 

End  of  program.  Press  ENTER  to  continue 

The  last  set  of  information  following  “Number  of  points  per  wave  length 
distribution  (%)”  is  the  resolution  data  for  this  geometry,  calculated  from  the  grid- 
network  and  wave  period  specified  by  the  user  in  the  general  information  file.  As  stated 
earlier,  less  than  5  grid  points  per  wavelength  is  not  valid,  over  15  points  per  wavelength 
is  generally  very  good.  The  resolution  data  shown  is  for  the  entire  computational  domain 
and  it  gives  the  resolution  distribution  for  all  the  elements.  This  information  is  output  to 
the  screen  only.  The  user  may  copy  and  retain  it  manually  for  later  reference.  This 
feature  is  also  useful  if  the  same  grid  is  to  be  used  again  with  a  different  wave  period;  the 
resolution  analysis  will  help  to  determine  whether  or  not  grid  modification  is  necessary. 
This  may  be  the  case  if  the  new  wave  period  is  smaller  than  the  original  design  period 
wave. 


60 


7.6  The  Wave  Model  -  Program  “cgwave” 

This  is  the  main  part  of  the  CGWAVE  model.  When  program  “cgwave”  is  run  for 
the  rectangular  harbor  case,  the  following  prompts  appear: 

Input  data  filename:  (*.dat) 
rect.dat 

Output  wave  potential  filename:  (*.out) 
rect.out 

A  block  of  information  reviewing  input  parameters  is  printed  to  the  screen.  This  is 
followed  by  messages  regarding  the  processing  of  input  information.  A  “hot  start”  (or  re¬ 
start)  option  then  follows: 

A  previous  solution  file  as  initial  state  (y/n)  ? 
n 

Iterations  Begin,  Please  Wait . 

The  model  will  run  until  either  the  convergence  parameter  or  the  maximum  number  of 
iterations  parameter  is  met  or  exceeded.  The  solution  will  then  be  written  to  the  file 
(*.out)  specified  earlier.  The  “hot  start”  option  allows  the  user  to  restart  the  model  with 
computations  beginning  from  an  old  CGWAVE  output  file.  This  can  be  useful  in  the 
event  that  the  model  is  forced  to  write  a  solution  (*.out)  before  the  convergence  criteria 
have  been  reached.  This  “partially  converged”  solution  will  occur  if  the  “Maximum 
Iterations  for  the  Linear  Equations”  parameter  is  exceeded  before  the  “Tolerance  for 
Convergence  of  Linear  Equations”  parameter  is  reached. 

7.7  Post-processing  the  Output  -  Program  “trans” 

The  output  from  program  “cgwave”  is  a  file  (e.g.  rect.out)  which  contains  grid 
information  and  calculated  complex  wave  elevations.  The  program  “trans”  may  be  used 
to  convert  these  potentials  into  a  readily  usable  format  or  to  obtain  output  at  desired 
locations. 
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7.7.1  Program  “trans”. 

Reformed  data  filename  (used  as  input  to  CGWAVE) :  (*.dat) 
rect.dat 

Wave  potential  file  (output  from  cgwave):  (*.out) 
rect.out 

Please  select  one  for  output: 
g  General  info  only 

t  General  info  with  complete  text  output 
n  None  above 

The  “General  info  only”  option  prints  text  information  about  the  input  parameters  and 
number  of  iterations  required  to  obtain  the  solution.  The  “General  info  with  complete 
text  output”  option  adds  the  node  numbers,  location,  computed  wave  amplitude  and  phase 
information  to  the  general  information  output.  The  naming  convention  used  for  the  file  is 
“filename.txt”. 

Create  output  of  amplitude  and  phase  for  SMS  (yes/no)  ? 

This  output  provides  information  in  a  format  that  is  readable  by  SMS  and  is  used  for 
graphic  viewing  of  amplitude,  phase  and  sea  surface  elevation  at  time  =  0.  The  naming 
convention  is  “filename.sol”.  This  convention  deviates  from  the  SMS  *.sol  file  naming 
convention  in  that  the  CGWAVE  *.sol  files  are  written  in  ASCII  text  format. 

Create  output  of  water  particle  kinematics  (yes/no)  ? 

This  output  provides  maximum  horizontal  velocities  and  pressure  gradients.  The  naming 
convention  is  “filename-kin.sol”. 

Create  time  series  output  file  for  velocity 
and  surface  elevation  animation  (yes/no)  ? 
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An  animated  view  of  the  velocity  and  surface  elevation  can  be  created  for  viewing  in 
SMS.  The  file  naming  convention  is  “filename.alp”. 

Do  you  want  output  at  selected  locations  (yes/no)  ? 

Amplitude  and  phase  data  can  be  extracted  at  discrete  output  points  or  along  lines 
(transects)  or  in  rectangular  sub-regions.  The  coordinates  of  the  points,  lines  and  sub- 
regions  may  be  specified  directly  from  the  keyboard  or  read  from  an  input  file.  The  input 
file  is  most  useful  for  collecting  data  at  many  points,  as  the  process  of  entering  all  the 
data  point  information  by  hand  can  become  tedious.  The  naming  convention  for  the 
output  is  “filename.sel”.  The  following  section  addresses  the  use  of  the  “output  at 
selected  locations”  option. 


7.7.2  Selected  Output  at  Points/Lines/Rectangular  Sub-regions 


As  stated  above,  there  are  two  available  methods  for  specifying  output  locations; 
by  keyboard  input  and  by  reading  from  a  file. 

Keyboard  input  Table  1  summarizes  the  information  needed  for  the  “keyboard  input” 
option  for  extracting  CGWAVE  predictions  at  points,  along  lines  and  in  rectangular  sub- 
regions.  Results  may  be  obtained  for  more  than  one  type  of  location  in  a  given  “trans” 
run. 


Table  1.  Summary  of  Input  Information,  Keyboard-Selected  Output  Option 


Method  of  Specifying  Location 

Keyboard  Input  Information  Needed  by 
Program  “trans” 

Specific  nodes 

Total  number  of  nodes,  node  numbers 

(x,y)  locations  (points) 

Total  number  of  points,  (x,y) 
coordinates 

Lines  (data  collected  along  transects) 

Total  number  of  lines,  (x,y)  location  of 
start  point,  (x,y)  location  of  end  point, 
number  of  intervals*  between  points. 
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Rectangular  sub-regions 


(x,y)  location  of  comer  points  in 
clockwise  order  beginning  from  upper 
left  point,  number  of  intervals*  between 
points  1  and  2  and  between  points  3  and 
4. 


♦Number  of  intervals  between  points  may  be  specified  as  different  values  for  each  line. 


Program  “trans”  will  ask  the  user  a  “yes/no”  question  regarding  each  method  of 
specifying  the  output  location.  A  “no”  answer  skips  to  the  next  method.  A  “yes”  answer 
will  result  in  additional  prompts  for  input  information.  After  the  input  is  complete,  the 
program  progresses  to  the  next  method. 


Input  File  Format.  The  same  output  location  data  given  in  Table  1  is  also  used  for  the 
“input  file”  option.  The  input  file  format  is  given  below.  The  naming  convention  for  this 
file  is  “filename.pls”.  The  “%”  symbols  are  comments  and  are  not  read  by  program 
“trans”,  however,  lines  beginning  with  “%”  should  not  split  the  individual  data  blocks. 
The  meaning  of  each  number  is  described  in  the  data  block.  The  sample  file  listed  below 
includes  all  available  output  options,  denoted  by  their  flags  “&id”  for  node  Ids,  “&pt”  for 
(x,y)  points,  “&ln”  for  transects,  “&rc”  for  rectangular  sub-regions.  The  option  blocks  do 
not  need  to  appear  in  the  order  shown,  nor  do  all  option  blocks  need  to  be  present.  For 
example,  if  the  user  wishes  to  specify  data  output  from  a  single  (x,y)  location,  then  the 
input  file  need  only  contain  the  “&pt”  flag  on  the  first  line,  the  quantity  of  points  (1  in 
this  example)  on  the  second  line,  the  x,y  location  of  the  point  on  the  third  line  and  the 
“end  of  data”  statement  on  the  last  line.  Note  that  the  “end  of  data”  statement  must  be  the 
last  line  of  the  file. 


%  Top  of  File 

% 

%  This  file  may  be  named  *.pls,  for  point,  line  and  square. 

% 

%  This  file  provides  information  for  selected  output  locations. 

% 
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%  Circular  island  problem 

% 

%  Selected  id  flag 
%  Total  number  of  node  ids 
%  id  numbers 
&id 
5 

10  23  78  55  36 
% 

%  Selected  points  flag 
%  Total  number  of  points 
%  x  y  location  of  each  point 
&pt 

3 

-8972.0  -6884.0 
3310.0  12353.0 
463.0  14456.0 

% 

%  Selected  Lines  flag 

%  Total  number  of  lines 

%  x  y  of  start  point  x  y  of  end  point 

%  number  of  intervals  between  points 

&ln 

3 

10000.0  0.0  37000.0  0.0 
29 

0.0  10000.0  0.0  37000.0 
29 

-10000.0  0.0  -37000.0  0.0 
29 
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% 

%  Rectangular  Sub-domain  flag 
%  upper  left  comer 
%  upper  right  comer 
%  lower  left  comer 
%  lower  right  comer 

%  number  of  intervals  between  first  two  comers 
%  number  of  intervals  between  last  two  comers 
&rc 

-4.0  4.0 
4.0  4.0 
4.0  -10.0 
-4.0  -10.0 
9 

29 

end  of  data 

7.8  Viewing  Solutions  in  SMS 

Any  of  the  *.sol  files  created  by  program  “trans”  contain  wave  information  that 
can  be  read,  viewed  and  printed  by  SMS.  To  do  so,  open  SMS  and  then  open  the 
corresponding  *.geo  file.  The  “Data  Browser”  tool  (also  found  under  the  “Data”  menu)  is 
used  to  import  the  solution  file.  After  opening  the  “Data  Set  Browser”  window,  select  the 
“Import”.  When  the  “Import  Data  Set”  window  appears,  select  “Generic  file”.  Choose 
the  appropriate  *.sol  file  from  the  “Open”  window.  When  viewing  Amplitude/Phase 
solutions,  three  time  step  values  will  appear  in  the  “Time”  window  (to  the  right  of  “Scalar 
data  set”)  in  the  data  set  browser.  The  time  value  “0.0”  corresponds  to  wave  amplitudes 
(which  have  values  greater  than  or  equal  to  0.0).  Time  “0.5”  corresponds  to  wave  phases 
(values  from  -1.0  to  1.0).  Time  “1.0”  corresponds  to  sea  surface  elevation  (both  positive 
and  negative  values).  After  choosing  the  type  of  solution  to  view,  click  the  “Done” 
button.  The  solution  contours,  or  vectors  if  vector  data  has  been  selected,  will  appear 
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within  the  modeled  domain.  Contours  and  vectors  can  be  manipulated  using  the  “Display 
Options”  tool  to  provide  a  meaningful  image. 
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7.9  Instructions  for  CGWAVE  Interface  in  SMS 
(MODE-2  Usage) 

7.9.1  Creating  a  Size  Function 

1.  Import  the  .xyz  data: 

File  Menu:  choose  Import 

Select  the  XYZ  data  button. 

Click  OK. 

Find  and  select  an  xyz  data  file  (bp.xyz,  for  example),  click  OPEN. 

2.  In  the  Mesh  Module  follow  these  steps  in  sequence: 

Data  Menu:  choose  Switch  Current  Model. 

Select  CGWAVE. 

Click  OK. 

CGWAVE  Menu:  choose  Create  Functions. 

Leave  only  the  Wavelength  option  on  (turn  all  others  off). 

Enter  desired  value  for  period. 

Click  OK. 

Data  Menu:  choose  Data  Calculator 

Select  a.  Wavelength  in  text  window  (upper  left). 

Click  Add  to  Expression. 

Select  “/”  (divide  symbol)  in  the  calculator. 

Enter  desired  number  of  elements  in  each  wavelength  (i.e.,10)  as 
denominator. 

Change  the  title  “new  data  set”  to  “size”  (i.e.  name  the  newly  created 
function). 

Click  on  Compute. 

Click  Done. 
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7.9.2  Creating  the  Background  Scatter  Set 

1 .  Create  Scattered  Data  Set 

Data  Menu:  choose  Mesh->Scattered  Data 

Click  OK  -  the  name  “scatter”  is  assigned  to  the  set  you  just  created. 

2.  Select  Interpolation  Method. 

Switch  to  Scatterpoint  Module. 

Interpolation  Menu:  choose  Interp  Opts. 

Choose  Inverse  Distance  Weighted. 

Click  Options. 

Choose  Constant  (Shepherd’s  Method). 

Click  OK  (Options  Dialog). 

Click  OK  (Interpolation  Method  Dialog). 

7.9.3  Creating  the  Coastline  and  Inlet  Arcs 

1 .  Import  the  Coastline  Arc. 

Go  to  the  Map  Module  and  choose  Feature  Object s\ Coverages. 

Change  the  Coverage  type  on  the  bottom  right  of  the  dialog  to  CGWAVE. 

Go  to  FileUmport  and  choose  the  Coastline  button. 

Find  the  coastline  file  and  push  OPEN. 

A  brown  Coastline  type  arc  will  be  created. 

2.  Create  the  Inlet  Arc. 

Zoom  in  on  the  harbor  area  of  the  arc  using  the  Zoom  tool. 

Use  the  Create  Feature  Arc  tool  to  create  an  arc.  The  arc  should  be  constructed 
such  that  it  spans  the  harbor  entrance  and  contains  a  single  vertex  located  near 
the  center  of  the  harbor  entrance  (the  arc  will  be  made  up  of  two  end  points  and 
a  vertex).  Note:  zoom  on  the  Inlet  part  of  the  harbor  and  make  sure  that  the  end 
nodes  of  the  new  arc  are  placed  exactly  on  one  of  the  existing  coastline  arc 
vertices  on  the  South  and  North  sides  of  the  inlet. 
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Using  the  Select  Feature  Arc  tool,  double  click  on  the  new  arc. 

Choose  Inlet  type  when  the  arc  attributes  dialog.  Press  OK  to  exit  the  dialog. 

Deselect  the  arc  by  clicking  on  a  blank  area  of  the  screen  (it  is  necessary  to 
deselect  any  selected  arcs  before  building  polygons). 

7.9.4  Defining  Domains 

1 .  Build  Polygons  and  Define  the  Harbor  Domain 

In  the  Feature  Objects  menu,  select  Build  Polygons. 

Use  the  Select  Feature  Polygons  tool  to  double  click  on  the  harbor  area  (which  is 
a  new  polygon)  and  change  the  CGWAVE  type  to  Harbor. 

2.  Defining  the  Domain  from  Coastline/Open  Boundary  Intersection  Points. 

Use  the  Select  Feature  Vertex  tool  to  select  a  vertex  on  coastline  corresponding  to 
the  point  where  the  coastline  and  open  boundary  will  meet.  If  no  vertex  exists 
at  this  location,  an  new  vertex  may  be  created  using  the  Create  Vertex  tool. 

With  the  vertex  selected,  go  to  Feature  Objects  menu:  choose  Vertices  <->Nodes. 

Repeat  step  2  for  the  remaining  coastline-open  boundary  intersection  point. 

Switch  to  Select  Feature  Point  /Node  tool. 

Click  the  start  point  of  the  open  boundary.  Hold  the  SHFT  key  and  click  the  end 
point  of  the  open  boundary.  Keep  in  mind  that  the  open  boundary  will  be 
created  in  a  counter-clockwise  direction. 

In  the  Feature  Objects  menu:  choose  Define  Domain. 

Enter  the  desired  Domain  Type  and  OK  to  exit. 

3.  Defining  the  Domain  from  the  Harbor. 

With  the  Select  Feature  Vertex  tool,  select  the  vertex  in  the  middle  of  the  Inlet 
(red)  arc. 

Go  to  Feature  Objects\Define  Domain  and  enter  the  desired  values.  Push  OK.  If 
a  prompt  appears  saying  that  it  was  unable  to  create  the  domain,  it  will  be 
necessary  to  change  the  settings  in  the  Define  Domain  dialog. 
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7.9.5  Assigning  Reflection  Coefficients 

Choose  Select  Feature  Vertex  tool. 

Select  vertices  on  the  coastline  where  changes  in  the  reflection  coefficients  are 
desired. 

Go  to  Feature  Objects\Vertics<->  Nodes.  This  splits  the  coastline  arc.  Each 
coastline  arc  can  be  assigned  its  own  reflection  coefficient.  (Multiple  vertices 
can  be  selected  and  converted  to  nodes  at  the  same  time  by  using  the  Shift  key 
to  multiply  select  the  vertices). 

Double  click  on  each  coastline  arc  and  specify  the  desired  reflection  coefficient 
for  that  segment.  (If  more  than  one  arc  will  have  the  same  reflection 
coefficient,  the  arcs  can  be  multiply  selected  using  the  Shift  key  and  the 
attributes  assigned  using  the  Feature  Objects\Attributes  command.) 

7.9.6  Assigning  Open  Ocean  Attributes 

Choose  Select  Feature  Arc  tool. 

Select  one  or  more  open  boundary  arcs. 

Select  Feature  Objectsi Attributes  command  or  double  click  on  a  single  arc  and 
specify  the  arcs  as  Open  Ocean  type.  (Note:  if  a  reflection  coefficient  is  desired 
along  an  ocean  boundary  to  simulate  a  near  coastline  type  situation,  the  ocean 
boundary  can  be  defined  as  a  coastline.  Use  the  instructions  for  assigning 
coastline  reflection  coefficients). 

7.9.7  Setting  the  Active  Scatter  Data  Function 

Choose  Feature  Objects\Build  Polygons. 

Double  click  on  the  domain  area  to  select  the  open  ocean  area  of  the  domain  and 
activate  the  attributes  dialog. 

Assign  it  to  be  an  Ocean  type. 

Select  the  Density  radio  button  and  click  on  the  Options  button. 

Select  the  “size”  function  (created  in  section  7.9.1)  for  Element  Size  and  the 
“elevation”  function  for  bathymetry. 
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Double  click  on  the  Harbor  polygon  and  repeat  Density  operation. 


7.9.8  Creating  a  Mesh  from  Feature  Objects 

De-select  all  polygons  by  clicking  somewhere  else  on  the  screen. 

Go  to  Feature  Objects\Map->2D  Mesh  and  push  OK  at  the  prompt. 

Toggle  the  Display  Meshing  Process  on  if  desired  and  push  OK. 

When  meshing  is  completed,  go  to  the  Scatter  Module  and  push  Display\Display 
Options  and  toggle  the  scatter  point  symbols  off  for  better  viewing  of  the  mesh 
(if  desired). 


7.9.9  Saving  the  CGWAVE  .dat  and  .xyz  (1-d)  files 

Go  to  the  Mesh  Module. 

Click  the  Select  Nodestrings  button  and  select  a  nodestring  (this  is  necessary  for 
renumbering). 

Choose  ElementslRenumber  and  push  OK.  Push  OK  if  prompted  again. 

Go  to  the  CGWAVE  menu,  choose  Save  Simulation. 

Select  file  names  for  the  .dat  and  .xyz  (1-d)  files.  If  it  is  indicated  that  the  1-d  file 
cannot  be  saved,  go  to  CGWAVE\Model  Control  and  increase  the  1-d  spacing 
or  number  of  1-d  nodes. 

7.9.10  Checking  Mesh  Quality 

Go  to  Scatter  Module. 

Choose  Data\Data  Browser  and  select  the  size  function  used  to  generate  the 
network. 

Choose  Inerpolate I  ..to  Mesh  and  enter  name  of  function  for  the  new  network  to 
represent  the  target  size  for  each  element. 

Go  to  Mesh  Module. 

Choose  CGWAVE\Create  Functions  and  turn  off  all  functions  except  grid  spacing. 

Click  OK. 
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Choose  Data\Data  Calculator  and  build  size  -  grid  function,  call  it  “error”.  This 
function  should  be  close  to  zero.  Areas  with  large  values  indicate  node  spacing 
that  does  not  match  the  target  spacing.  If  the  spacing  is  too  far  off,  contact 
BYU  or  edit  mesh  by  hand.  Another  function  useful  for  assessing  mesh  quality 
is  the  error/size,  call  it  “pcnt_err”.  Areas  with  red  color  point  to  the  elements, 
which  are  bigger  than  the  target  size,  so  refinement  may  be  warranted. 

7.9.11  Running  CGWAVE 

CGWAVE  can  be  run  from  CGWAVElRun  CGWAVE.  Click  on  the  file  button  to 
give  the  directory  that  contains  the  CGWAVE  executable  file. 


73 


8  EXAMPLES 


Surface  water  waves  have  been  modeled  using  CGWAVE  for  a  number  of 
idealized  testcases  and  real  practical  applications.  A  number  of  these  applications  are 
presented  here  to  both  validate  CGWAVE  and  to  demonstrate  its  ability  to  model  waves. 

8.1  Shoal  on  a  Sloping  Beach 

Surface  water  waves  propagating  over  the  shoal  presented  by  Berkhoff,  et  al. 
(1982)  were  modeled.  This  test  case  demonstrates  the  ability  of  CGWAVE  to  simulate 
the  effects  of  complex  coastal  bathymetric  features.  The  shoal  in  Figure  7  is  oriented 
such  that  the  major  axis  of  the  shoal  is  parallel  to  contours  of  the  water  depth.  Data  is 
collected  along  the  eight  sections  (transects)  shown  in  Figure  7. 

Normal  incident,  plane  waves  having  a  period  of  1  second  are  modeled  for  all 
runs.  An  incident  amplitude  of  1.0  m  is  used  for  model  runs  using  the  linear  dispersion 
relation.  Runs  that  are  made  using  the  non-linear  dispersion  relation  use  an  incident 
amplitude  of  0.0232  m. 

Results  were  compared  for  model  runs  made  with  different  grid  resolution.  This 
was  accomplished  by  vaiying  the  number  of  nodes  in  the  computational  domain  and  thus 
constructing  different  finite  element  grids.  Nodes  were  located  such  that  the  number  of 
nodes  per  wavelength  remains  constant  throughout  the  domain.  Grid  densities  from  three 
to  fifteen  elements  per  (local)  wavelength  were  used.  The  resulting  finite  element  grids 
contain  between  2500  and  75000  nodes  and  5000  and  150000  elements,  respectively. 
CPU  time  is  not  taken  into  account  when  comparing  the  quality  of  the  model  results,  as 
the  longest  solution  run  is  completed  in  less  than  four  hours  on  a  desktop  (200  MHz 
processor)  PC. 
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Figures  8,  9  and  10  show  results  for  domains  having  resolutions  of  3,  6,  8  and  15 
nodes  per  wavelength.  The  numerical  solution  in  Figure  8  shows  little  correlation  to  the 
lab  data.  The  results  in  Figures  9  and  10  show  a  progression  to  the  lab  data  as  the 
resolution  increases.  The  linear  solution  in  Figure  10  matches  the  data  better  than  either 
of  the  previous  resolutions,  however  section  7  does  not  agree  with  the  data  in  the  region  6 
to  10  meters  behind  the  shoal.  A  grid  resolution  of  about  10  or  greater  nodes  per 
wavelength  is  generally  thought  to  be  necessary  for  an  accurate  solution  to  most 
problems. 

The  non-linear  dispersion  mechanism  was  applied  only  to  grids  containing  15 
nodes  per  wavelength.  This  mechanism  was  not  run  for  lower  resolution  grids  because 
the  linear  amplitude  solution  is  used  in  the  formulation  of  the  non-linear  effect.  Inclusion 
of  the  non-linear  dispersion  relation  in  the  calculations  has  a  noticeable  effect  on  the 
solution  for  all  sections  (Figure  10).  Sections  3  and  5  show  a  closer  fit  to  the  data  in  the 
areas  to  the  left  and  right  of  center.  Section  7  shows  the  most  marked  improvement  in  the 
region  that  begins  6  to  7  m  behind  the  shoal. 

This  example  problem  showed  that  the  results  produced  by  the  linear  run  of  CGWAVE 
are  good,  provided  that  a  fine  enough  grid  resolution  is  used.  The  best  results  are 
obtained  with  the  highest  grid  resolution.  The  addition  of  the  non-linear  dispersion 
relation  (wave-wave  interactions)  mode  shows  a  substantial  improvement  in  the  model 
estimates,  yielding  the  best  comparison  to  the  measurements. 

59-83. 
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Figure  7.  Modeled  domain  of  Berkhoff  et  al.  (1982) 
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Figure  9.  Wave  height  comparisons  for  increasing  domain  resolutions 
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Figure  10.  Wave  height  comparisons  from  linear  and  non-linear  CGWAVE  runs  on  a 
domain  resolution  of  15  points  per  wavelength 
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8.2  Irregular  Wave  Propagation  Over  a  Shoal 

The  hydraulic  model  study  of  Vincent  and  Briggs  (1989)  was  used  to  simulate 
spectral  wave  conditions  in  CGWAVE.  The  directional  spectral  wave  generator  (DSWG) 
at  the  Coastal  Engineering  Research  Center  of  the  U.S.  Army  Engineers  was  used  to 
simulate  refraction-diffraction  of  directionally  spread  irregular  waves  over  an  elliptic 
shoal.  This  study  showed  that  spectral  conditions  can  produce  significantly  different 
results  than  monochromatic  waves.  The  model  domain  consists  of  an  elliptic  shoal 
surrounded  by  a  region  of  constant  depth.  The  orientation  of  the  shoal  and  placement  of 
coastal  boundaries  are  shown  in  Figure  11.  The  minimum  water  depth  over  the  shoal  is 
approximately  0.15  m.  The  depth  around  the  shoal  is  a  constant  0.4572  m.  Wave  height 
data  was  collected  along  the  nine  sections  indicated  in  Figure  11. 

The  spectral  input  to  CGWAVE  was  constructed  using  a  numerically  generated, 
two-dimensional  spectra,  E(f,0)  =  E(f)D(0).  Expressions  used  for  the  energy  spectrum, 
E(f),  and  the  directional  spectrum,  D(0),  follow.  The  energy  spectrum  is  obtained  using 
the  Texel  Marsen  Arsloe  (TMA)  spectrum  (Bouws  et  al.  1985,  Panchang  et  al.  1990), 

where  a  =  Phillips’  constant 
fm  =  peak  frequency 

y  =  peak  enhancement  factor  (=1  for  the  Pierson-Moskowitz  spectrum) 

G  =  shape  parameter  (G  =  Ga  if  f  <  fm and  G  =  Gb  if  f  >  fm) 


'f.'4 


E (/)  =  g2a(27r)~  /'5  exp<|  -  1.25|^J  +  (in  r)exp| 


jf-fj 

2  o2f2m 


<| )(f,h)  is  a  function  that  incorporates  depth  effects  and  is  approximated  by  (Hughes  1984) 


as 


f 0.5(g))2  for  coh  <  1 


((,: 


1  —  0.5(2  —  CDh)2 

1.0 


for  1  <  coh  <2 
for  (Dh  >  2 


where  C0h  =  2ttf(h/g)1/2. 
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The  directional  spectrum  is  calculated  as  follows  (Vincent  and  Briggs  1989), 

(M,)2 


D(0)  =  ^I  +  7Sex  P 

27T  K  j=1 


cos  j(e-9n) 


where  0m  =  mean  wave  direction  =  0° 

j  =  number  of  terms  in  the  series  (20  in  numerical  calculations) 
cm  =  spreading  parameter 


The  same  functions  E(f)  and  D(0)  were  also  used  in  the  laboratory  study  by 
Vincent  and  Briggs  (1989)  to  construct  a  number  of  frequency  and  directional  spectra  (for 
details,  see  Vincent  and  Briggs  1989  or  Panchang  et  al.  1990).  One  spectral  case,  Bl, 
consisted  of  a  broad  frequency  spectrum;  the  other  case,  Nl,  had  a  narrow  frequency 
spectrum.  The  input  data  are  summarized  in  Table  2.  Two  monochromatic  and  two 
spectral  cases  are  shown.  The  incident  wave  height  is  2[2E(/)D(0)A/A0] 1/2 . 


Table  2.  Model  Input  Data 


Input 

Case  ID 

Period 

(sec) 

Significant 

Wave 

Height  (m) 

a 

■ 

Monochromatic 

Ml 

1.3 

.0550 

— 

P  ■ 

M2 

1.3 

.0254 

~ 

— 

— 

Broad-directional 

Bl 

1.3 

.0775 

.01440 

2 

30° 

Narrow-directional 

Nl 

1.3 

.0775 

.01440 

2 

10° 

The  discretization  of  the  directional  spectrum  is  summarized  in  Table  3.  Other 
researchers  have  used  the  following  discretizations.  Panchang  et  al.  (1990)  used  a 
directional  bin  width  (A0)  of  4.39°  in  a  range  from  45°  to  +45°,  resulting  in  39  directional 
components  for  all  spectra.  Zhao  and  Anastasiou  (1993)  found  that  as  many  as  63 
components  were  needed  for  the  Bl  and  Nl  cases.  Li  et  al.  (1993)  produced  results  that 
suggest  that  an  adequate  solution  can  be  obtained  with  as  few  as  five  directional 
components  and  five  frequency  components.  In  light  of  the  results  of  Li  et  al.  (1993),  it 
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was  decided  that  some  exploration  of  the  necessary  number  of  directional  components 
was  in  order. 


Table  3.  Directional  Resolutions  Used  for  Model  Input 


Case  ID 

Number  of  Directional 
Components  (bins) 

A0(°) 

B1 

9 

20 

29 

4.14 

N1 

5 

12 

15 

4 

Tests  were  also  performed  to  quantify  the  effects  of  the  boundary  conditions  on 
the  solution.  It  was  hypothesized  that  reflections  from  the  coastal  boundaries  (simulated 
walls  of  the  wave  tank)  and  improper  treatment  of  these  fully  absorbing  boundaries  could 
cause  significant  variations  in  the  final  solution,  especially  for  components  with  large 
angles  of  incidence.  This  hypothesis  was  tested  by  comparing  the  results  obtained  from 
two  domains.  Domain  A  closely  approximated  the  physical  model  and  contained  coastal 
boundaries  (Figure  11).  Interior  and  exterior  boundaries  were  considered  fully  absorbing 
(Kr=  0.0  and  KExt  =  0.0).  Domain  B  used  a  full  circle  as  an  open  boundary  and 
contained  all  the  depth  variations  (Figure  12).  It  was  found  that  the  boundary  conditions 
do  not  significantly  effect  the  results  from  CGWAVE  in  the  area  of  interest.  Results  are 
only  presented  for  Domain  A. 

CGWAVE  results  for  the  monochromatic  conditions  Ml  and  M2  are  shown  in 
Figures  13  and  14.  The  wave  amplitudes  double  behind  the  shoal  and  decrease  on  the 
sides  of  the  shoal.  The  linear  model  results  show  a  similar  pattern.  The  non-linear 
numerical  calculations  match  the  data  better  than  the  purely  linear  results.  The  data 
shows  that  the  maximum  wave  height  occurs  further  behind  the  shoal  than  predicted  by 
both  linear  and  non-linear  solutions  (section  7,  Figures  13  and  14). 

The  CGWAVE  results  for  the  B1  spectral  wave  condition  with  A0  =  20° 
discretization  is  shown  in  Figure  15.  The  result  exhibits  finger-like  areas  of  large  and 
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small  wave  height  radiating  away  from  the  shoal  in  the  down-wave  direction  that  are  not 
present  in  the  lab  data.  The  model  results  in  Figure  15  do  not  show  the  same  level  of 
modulation  as  the  lab  results.  The  N1  solutions,  although  not  shown,  exhibit  the  same 
down  wave  characteristics  as  the  B1  solutions. 

The  numerical  solution  approaches  the  lab  data  when  the  directional  resolution  is 
increased.  Figures  16  and  17  show  the  results  for  the  B1  and  N1  spectra,  respectively, 
using  A0  equal  to  4.14°.  The  widths  of  the  directional  spectra  were  limited  to  +/-  60°  for 
B1  and  +/-300  for  N1  spectra,  since  negligible  amounts  of  energy  are  outside  this  range. 
The  finger  like  patterns  found  in  the  cases  with  large  A0  were  not  seen  in  these  results. 
This  finding  agrees  with  the  observations  made  by  Zhao  and  Anastasiou  (1993)  and 
Panchang  et  al.  (1990)  regarding  the  width  of  directional  bins  needed  in  the  input  spectra. 

The  effect  of  amplitude-dependent  non-linear  dispersion  was  examined  for  the  B 1 
and  N1  cases.  Zhao  and  Anastasiou  (1993)show  an  improvement  in  the  numerical 
solution.  The  CGWAVE  results  are  essentially  unchanged  (Figures  16  and  17).  This 
occurs  because  CGWAVE  operates  on  each  spectral  component  independently.  Once  a 
linear  solution  for  a  given  component  is  found,  these  solution  amplitudes  are  used  in  the 
non-linear  dispersion  relation  and  the  new  solution  is  found.  The  new  amplitudes  are 
again  used  in  the  non-linear  mechanism.  These  iterations  continue  until  a  convergence 
criterion  is  reached.  The  final  solution  to  the  component  is  saved  and  the  model  is  run  for 
the  next  component.  The  incident  wave  heights  associated  with  each  spectral  component 
are,  at  their  largest,  an  order  of  magnitude  smaller  than  the  significant  wave  height.  The 
wave  heights  used  to  compute  non-linear  effects  are  only  a  fraction  of  the  significant 
wave  height,  causing  the  non-linear  dispersion  relation  to  have  a  negligible  effect. 

For  this  problem,  CGWAVE  produces  linear  results  that  closely  match  the  data 
for  the  monochromatic  waves  propagating  over  the  shoal.  The  inclusion  of  wave- wave 
interaction  has  improved  the  quality  of  this  solution.  CGWAVE  also  produces  linear 
results  that  closely  match  the  data  for  the  irregular  wave  cases  when  adequate  directional 
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resolution  is  provided  in  the  input  spectra.  The  solution  algorithm  used  in  the  current 
version  of  the  model  does  not  completely  allow  non-linear  effects  to  be  accurately 
represented  for  the  spectral  wave  modeling.  Other  solution  algorithms  should  be 
evaluated  in  order  to  determine  whether  or  not  they  may  yield  results  that  more  closely 
match  the  lab  data. 
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Figure  1 1 .  Model  domain  for  wave  propagation  over  an  elliptic  shoal  [after  Vincent  and 
Briggs  (1989)] 


Figure  12.  Domain  without  coastal  boundary  effects 


Figure  13.  Wave  height  comparisons  for  Ml  input  condition 
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Figure  14.  Wave  height  comparisons  for  M2  input  condition 
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Figure  15.  Estimated  wave  heights  for  B1  input  spectrum,  20°  directional  discretization 
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Figure  16.  Wave  height  comparisons  for  B1  input  condition,  4.14°  directional 
discretization 
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Figure  17.  Wave  height  comparisons  for  N1  input  condition,  4.14°  directional 
discretization 
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8.3  Rectangular  Harbor  -  Resonance 

In  this  example,  we  examine  the  influence  of  bottom  friction  on  dispersion  of 
wave  energy  by  friction  for  resonating  waves  in  a  rectangular  harbor.  The  open  boundary 
and  harbor  dimensions  are  shown  in  Figure  18.  Theoretical  and  lab  data  for  this  case 
were  presented  in  Ippen  and  Goda  (1963),  Lee  (1971)  and  Chen  (1986).  Ippen  and  Goda 
(1963)  showed  that  linear  theory  predicts  amplitudes  that  are  too  large  near  the  resonant 
frequency  of  the  harbor.  Chen  (1986)  examined  the  effects  of  bottom  friction  and  coastal 
reflection  on  the  amplitudes.  Laboratory  data  was  collected  at  the  center  of  the  back  wall 
of  the  harbor  (Lee,  1969). 

Model  input  conditions  were  obtained  from  Lee  (1969)  and  Chen  (1986).  Table  4 
summarizes  the  input  wave  conditions. 

Table  4.  Model  Input  Data 


Amplitude 

k  l 

Amplitude 

K  / 

.13 

0.2 

.30 

1.8 

.13 

0.4 

.25 

2.0 

.13 

0.6 

.25 

2.5 

.13 

0.8 

.30 

3.0 

.13 

1.0 

.25 

3.5 

.13 

1.1 

.35 

4.0 

.38 

1.2 

.35 

4.2 

.38 

1.3 

.35 

4.4 

.25 

1.4 

.25 

4.6 

.25 

1.5 

.30 

4.8 

.25 

1.6 

.30 

5.0 

The  term  k  is  the  wave  number  and  /  (0.31 1 1  m)  is  the  length  of  the  harbor.  Variations  in 
the  friction  coefficient  and  coastal  reflection  were  considered.  Waves  were  normally 
incident  to  the  exterior  boundary  for  all  model  runs. 

Results  for  a  fully  reflecting  harbor  are  shown  in  Figures  19,  20  and  21.  The 
results  of  linear  CGWAVE  run  without  friction  show  a  good  match  to  the  analytical 
solution  of  Ippen  and  Goda  (1963)  (Figure  19).  Figure  20  shows  that  the  amplification  at 
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the  resonant  peaks  decreases  as  the  coastal  boundary  reflection  coefficient  decreases. 
These  reflection  coefficient  does  not  affect  the  frequencies  at  which  the  resonant  peak 
occurs. 


Figure  21  shows  the  effects  of  friction  on  the  CGWAVE  predictions.  As  the 
friction  parameter  is  increased,  the  solution  shows  a  progression  toward  the  lab  data  at  the 
first  resonant  peak.  The  peak  frequency  shows  a  small  shift  for  the  larger  values  of  the 
friction  parameter.  At  the  second  resonant  peak  the  solutions  with  friction  match  the 
solution  without  friction,  indicating  that  energy  dissipation  due  to  bottom  friction  is  less 
important  for  short  period  waves  than  for  long  period  waves. 

For  this  problem,  CGWAVE  reproduces  laboratory  data  by  including  the  effects 
of  frictional  dissipation  and  boundary  absorption.  The  CGWAVE  results  show  that 
harbor  resonance  is  sensitive  to  the  energy  absorbed  at  the  coastal  boundaries  and  to 
dissipation  by  friction.  The  choice  of  friction  parameter  and  coastal  reflection  coefficient 
does  clearly  affect  the  model  estimates  in  resonance  cases. 
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Figure  18.  Rectangular  resonance  harbor  with  semi-circular  open  boundary 
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Figure  19.  Theoretical  and  numerical  resonance  curves  for  a  fully  reflecting  rectangular 
harbor 
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Figure  20.  Harbor  response  curves  for  various  values  of  coastal  reflection  coefficient 


Re  so  on  se  Curves  for  Rectangular  Harbor 


Figure  21.  Harbor  response  curves  for  various  values  of  the  friction  factor 
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8.4  Circular  Island 


In  this  example  problem,  we  consider  the  propagation  of  long  waves  around  a 
circular  island  is  simulated.  Homma  (1950)  and  Jonsson,  et  al.  (1976)  showed  that 
extensive  scattering  occurs  for  long  waves.  This  case  tests  and  demonstrates  the  ability  of 
CGWAVE  to  handle  strong  scattering  effects  with  a  relaxed  (Xu  et  al.  1996)  boundary 
condition.  This  boundary  condition  allows  the  model  to  be  applied  to  a  greater  range  of 
modeling  problems,  as  it  does  not  require  full  reflection  of  the  coastline  outside  the  model 
domain. 

The  modeled  domain  (Figure  22)  consists  of  a  circular  island  on  a  shoal  with  a 
parabolic  variation  in  depth.  The  island  has  a  radius  of  10  km.  The  open  boundary  is 
placed  at  a  radius  of  35  km  from  the  center  of  the  island.  The  minimum  water  depth  is 
444  m  at  the  edge  of  the  island  and  increases  to  4000  m  at  the  outer  edge  of  the  shoal. 
The  domain  has  constant  depth  beyond  the  shoal.  Incident  wave  periods  of  240,  410  and 
480  seconds  are  considered. 

Figure  23  shows  the  results  for  the  240  second  incident  wave  condition.  This 
figure  has  been  divided  in  half  due  to  symmetry  of  the  solution.  Waves  are  incident  from 
the  left.  The  upper  portion  is  the  analytical  solution,  the  lower  portion  is  the  CGWAVE 
result.  The  CGWAVE  results  match  the  analytical  solution.  Notably  absent  are  signs  of 
artificial  boundary  effects  near  the  open  boundary.  Similar  results  were  obtained  for  the 
longer  wave  periods,  with  results  matching  the  analytical  data. 
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INCIDENT  WAVE 


8.5  Ponce  de  Leon  Inlet,  FL 


This  example  shows  the  application  of  CGWAVE  to  Ponce  de  Leon  Inlet,  Florida. 
This  inlet  is  the  subject  of  a  physical  model  study,  allowing  the  comparison  of  laboratory 
data  to  CGWAVE  model  runs.  Wave-current  interaction  may  play  a  very  important  role 
on  the  wave  fields  to  be  predicted  in  coastal  inlets  (Demirbilek  et  al.  1996b);  this 
capability  will  be  available  in  CGWAVE  model  in  the  near  future. 

The  coastline  of  the  modeled  domain  is  oriented  in  roughly  a  north/south  manner. 
The  modeled  domain  is  bounded  by  a  semi-circle  of  2.5  km  radius.  A  jetty  extends  into 
the  domain  from  the  north  side  of  the  inlet.  The  finite  element  mesh  for  this  domain 
contains  of  about  118000  nodes  and  234000  elements. 

An  example  of  a  CGWAVE  run  of  this  domain  is  given.  The  wave  phase  diagram 
is  shown  in  Figure  24  and  estimated  wave  amplitudes  are  shown  in  Figure  25.  The 
incident  wave  amplitude  was  15  seconds.  The  incident  wave  direction  was  normal  to  the 
coastline.  Coastlines  were  chosen  to  be  fully  energy  absorbing  (reflection  coefficient  of 
0.0).  Results  from  the  numerical  model  are  compared  to  laboratory  data  in  Figure  26. 
The  data  comparison  is  made  along  a  transect  oriented  parallel  to  the  coastline  and 
located  several  hundred  meters  to  the  south  of  the  jetty.  The  CGWAVE  model  results  are 
for  this  case  are  shown  prior  to  fine  tuning  of  the  numerical  code  and  adjustment  of 
model  parameters.  The  study  is  in  progress  at  the  time  of  this  writing.  A  detailed 
discussion  of  the  results  is  expected  to  appear  in  a  relevant  scholarly  journal  at  the 
conclusion  of  the  study. 
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Figure  24.  Phase  diagram  for  normal  incident  direction,  15-second  period  wave  condition 
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2.50 


Figure  25.  Estimated  wave  height  diagram  for  normal  incident  direction,  15-second 
period  wave  condition 


Gauge  Locations  Along  Transect  (m) 


Figure  26.  Comparison  between  laboratory  data  and  CGWAVE  wave  height  estimates 
along  a  transect 
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8.6  Barbers  Point  Harbor,  HI 


This  example  shows  the  application  of  CGWAVE  to  the  harbor  and  near-by 
coastal  regions  of  Barbers  Point  Harbor,  Hawaii.  A  study  is  being  conducted  to  explore 
harbor  response  to  wave  periods  between  30  and  4000  seconds.  Of  particular  interest  are 
the  contributions  infragravity  waves  to  the  resonance  of  the  harbor. 

The  modeled  domain  consists  of  the  harbor  area  and  a  region  outside  the  harbor  of 
3  km  radius.  The  radius  of  the  area  outside  the  harbor  is  measured  from  the  mouth  of  the 
harbor.  This  domain  size  allows  multiple  wavelengths  to  be  modeled  within  the  domain 
for  all  but  the  longest  period  waves.  With  appropriate  resolution  (based  on  the  shortest 
period  wave)  the  grid  contains  33000  nodes  and  65000  elements.  Both  monochromatic 
and  spectral  conditions  are  used  as  input  in  this  modeling  study. 

An  example  of  the  estimated  wave  phases  and  amplitudes  is  given  in  Figures  27 
and  28,  respectively.  These  estimates  are  based  on  a  normal  incident  input  wave  with  a 
period  of  100  s  and  incident  amplitude  of  1  m.  The  coastline  has  been  assumed  to  be 
fully  absorbing  for  this  example  case.  A  detailed  analysis  of  the  results  of  this  study  is  in 
progress  at  the  time  of  this  writing.  Comparisons  of  CGWAVE  wave  height  estimates 
are  to  be  made  with  field  measurements  in  and  around  the  harbor  region.  The  results  and 
analysis  are  to  be  documented  in  relevant  journal  articles  upon  completion  of  the  study. 
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8.7  Onslo  Bay,  NC 


This  example  shows  the  application  of  CGWAVE  to  New  River  Inlet  in  Onslo 
Bay,  North  Carolina.  CGWAVE  was  applied  to  this  region  to  provide  wave  height 
estimates  for  JLOTS  (Joint  Logistics  Over-The-Shore)  exercises. 


A  rectangular  open  boundary  is  used  to  separate  the  modeled  domain  from  the 
open  ocean.  The  modeled  area  outside  the  inlet  is  3  km  along  shore  by  2  km  cross  shore. 
The  inlet  area  of  the  domain  extends  approximately  1.5  km  inland.  The  finite  element 
mesh  for  this  domain  contains  of  about  30000  nodes  and  60000  elements.  The  grid 
resolution  is  good  (no  less  than  10  nodes  per  wavelength)  for  waves  of  period  equal  to  or 
longer  than  35  seconds  are. 


The  phase  and  amplitude  results  for  a  monochromatic  wave  simulation  are  shown  in 
Figures  29  and  30,  respectively.  For  this  model  run,  incident  waves  were  normal  to  the 
coastline.  The  incident  wave  period  was  35  seconds.  The  incident  wave  height  was  1  m. 
The  coastline  was  assigned  a  reflection  coefficient  of  0.4.  The  alternating  pattern  of  wave 
phase  is  shown  in  Figure  29.  Wave  heights  of  2  m  were  predicted  to  the  north  side  of  the 
inlet,  represented  by  the  darkest  color  shown  in  Figure  30. 
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Figure  29.  Wave  phase  diagram  for  normal  incident  direction,  35-second  period  wave 
condition 


107 


Figure  30.  Estimated  wave  height  diagram  for  normal  incident  direction,  35-second 
period  wave  condition 
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